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ABSTRACT 

We show that the deterministic past history of the Universe can be uniquely recon- 
structed from the knowledge of the present mass density field, the latter being inferred 
from the 3D distribution of luminous matter, assumed to be tracing the distribution 
of dark matter up to a known bias. Reconstruction ceases to be unique below those 
scales - a few Mpc - where multi-streaming becomes significant. Above 6h~^ Mpc we 
propose and implement an effective Monge-Ampere-Kantorovich method of unique 
reconstruction. At such scales the Zel'dovich approximation is well satisfied and re- 
construction becomes an i nstance of optimal mass transportation, a problem which 
goes back tojMongc (17811) . After discretization into N point masses one obtains an 
assignment problem that can be handled by effective algorithms with not more than 
0{N^) time complexity and reasonable CPU time requirements. Testing against TV- 
body cosmological simulations gives over 60% of exactly reconstructed points. 

We apply several interrelated tools from optimization theory that were not used 
in cosmological reconstruction before, such as the Monge-Ampere equation, its re- 
lation to the mass transportation problem, the Kantorovich duality and the auction 
algorithm for optimal assignment. Self-contained discussion of relevant notions and 
techniques is provided. 

Key words: cosmology: theory - large-scale structure of the Universe - hydrody- 
namics 



1 INTRODUCTION 

Can one follow back in time to initial locations the highly 
structured present distribution of mass in the Universe, as 
mapped by redshift catalogues of galaxies? At first this 
seems an ill-posed problem since little is known about the 
peculiar velocities of galaxies, so that equations governing 
the dynamics cannot just be integrated back in time. In fact, 
it is precisely one of the goals of reconstruction to deter mine 
the pe culiar velocities. Since the pioneering work of Pe eblgsj 
|l989?), a number of reconstruction techniques have been 
proposed, which frequently provided non-unique answers.^ 
Cosmological reconstruction should however take ad- 
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^ We put the present work in context of several important exist- 
ing techniques in Section 171 



vantage of our knowledge that the initial mass distribu- 
tion was quasi-uniform at bary on-photon de coupling, about 
14 billion years ago (see, e.g.. [Sus^ereegi Binncv 1994|). 
In a recent Letter to Nature l)FYiscI^^ani2002f) . four of us 
have shown that, with suitable assumptions, this a priori 
knowledge of the initial density field makes reconstruction 
a well-posed instance of what is called the optimal mass 
transportation problem. 

A well-known fact is that, in an expanding universe with 
self-gravitating matter, the initial velocity field is 'slaved' to 
the initial gravitational field, which is potential; both fields 
thus depend on a single scalar function. Hence the number 
of unknowns matches the number of constraints, namely the 
single density function characterising the present distribu- 
tion of mass. 

This observation alone, of course, does not ensure 
uniqueness of the reconstruction. For this, two restrictions 



© 0000 RAS 



2 Y. Brenier et al. 




deblais remblais 

Figure 1. A sketch of Monge's mass transportation problem in 
which one searches the optimal way of transporting earth from 
cuts {deblais) to fills (remblais), each of prescribed shape; the 
cost of transporting a molecule of earth is a given function of the 
distance. The MAK method of reconstructing the early Universe 
described in this paper corresponds to a quadratic cost. 

will turn out to be crucial. First, from standard redshift 
catalogues it is impossible to resolve individual streams of 
matter with different velocities if they occupy the same space 
volume. This 'multi-streaming' is typically confined to rela- 
tively small scales of a few megaparsecs (Mpc), below which 
reconstruction is hardly feasible. Second, to reconstruct a 
given finite patch of the present Universe, we need to know 
its initial shape at least approximately. 

It is our purpose in the present paper to clarify the phys- 
ical nature of the factors permitting a unique reconstruction 
and of obstacles limiting it, and to give a detailed account 
of the way some recent developments in the optimal mass 
transportation theory are applicable. (Fig. may give the 
reader some feeling of what mass transportation is about.) 

The paper is organized as follows. In Section |21 we for- 
mulate the reconstruction problem in an expanding universe 
and state the main result about uniqueness of the solution. 

In the next three sections we devise and test a re- 
construction technique called MAK (for Monge-Ampere- 
Kantorovich) within a restricted framework where the La- 
grangian map from initial to present mass locations is taken 
potential. In Section |31 we discuss the validity of the poten- 
tiality assumption and its relation to various approximations 
used in cosmology; then we derive the Monge-Ampere equa- 
tion, a simple consequence of mass conservation, introduce 
its modern reformulation as a Monge-Kantorovich problem 
of optimal mass transportation and finally discuss different 
limitations on uniqueness of the reconstruction. In Section^] 
we show how discretization turns optimization into an in- 
stance of the standard assignment problem; we then present 
effective algorithms for its solution, foremost the 'auction' 
algorithm of D. Bertsekas. Section |3 is devoted to testing 
the MAK reconstruction against TV-body cosmological sim- 
ulations. 

In Sectional we show how the general case, without the 
potentiality assumption, can also be recast as an optimiza- 
tion problem with a unique solution and indicate a possible 
numerical strategy for such reconstruction. In Section Q we 
compare our reconstruction method with other approaches 
in the literature. In Section |H1 we discuss perspectives and 
open problems. 

A number of topics are left for appendices. In Ap- 
pendix ^ we derive the Eulerian and Lagrangian equa- 
tions in the form used throughout the paper (and provide 
some background for non-cosmologists) . Appendix^ is de- 
voted to the history of optimal mass transportation the- 



ory, a subject more than two centuries old llMongelll78l . 
which has undergone significant progress within the last two 
decades. Appendix El is a brief elementary introduction to 
the technique of duality in optimization, which we use sev- 
eral times throughout the paper. Appendix ^ gives details 
of the uniqueness proof that is only outlined in Sectional 

Finally, a word about notation (see also Appendix EJ- 
We are using comoving coordinates denoted by a; in a frame 
following expansion of the Universe. Our time variable is not 
the cosmic time but the so-called linear growth factor, here 
denoted by r, whose use gives to certain equations the same 
form as for compressible fiuid dynamics in a non-expanding 
medium. The subscript refers to the present time (redshift 
z = 0), while the quantities evaluated at the initial epoch 
take the subscript or superscript 'in.' Following cosmological 
usage, the Lagrangian coordinate is denoted q. 



2 RECONSTRUCTION IN AN EXPANDING 
UNIVERSE 

The most widely accepted explanation of the large-scale 
structure seen in galaxy surveys is that it results from small 
primordial fiuctuations that grew under gravitational self- 
interaction of coUisionless cold dar k matter (CDM) particles 
in an expanding universe (see, e.g.. lBernardeati et alJ (|2002h 
and references therein). The relevant equations of motion, 
derived in Appendix El sire the Euler-Poisson equations'^ 
written here for a flat, matter-dominated Ei nstein-de Sit- 
ter u niverse (for more general case see, e.g., ICatelan et alJ 

3 

drV + (V ■ V j:)v = - — (-U + Va;l^g), (1) 
AT 

drp + V, • {pv) = 0, (2) 
V^^g = (3) 

T 

Here v denotes the velocity, p denotes the density (normal- 
ized by the background density g) and i^g is a rescaled grav- 
itational potential. All quantities are expressed in comoving 
spatial coordinates x and linear growth factor r, which is 
used as the time variable; in particular, v is the Lagrangian 
r-time derivative of the comoving coordinate of a fluid ele- 
ment. 

2.1 Slaving in early-time dynamics and its fossils 

The right-hand sides of the momentum and Poisson equa- 
tions Q and ^ contain denominators proportional to r. 
Hence, a necessary condition for the problem not to be sin- 
gular as r — > is 

vin{x) + Vi^t^g" = 0, pm(a;) = 1. (4) 

In other words, (i) the initial velocity must be equal to (mi- 
nus) the gradient of the initial gravitational potential and 
(ii) the initial normalized mass distribution is uniform. We 
shall refer to these conditions as slaving. Note that the den- 
sity contrast p— 1 vanishes initially, but the rescaled gravita- 
tional potential and the velocity, as defined here, stay finite 

^ Also often called the Euler equations. 
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thanks to our choice of the hnear growth factor as time vari- 
able. Therefore we refer to the initial mass distribution as 
'quasi-uniform.' 

In the sequel, when we mention the Euler-Poisson 
initial-value problem, it is always understood that we start at 
T — and assume slaving. Hence we are extending the New- 
tonian matter-dominated post-decoupling description back 
to r = 0. By examination of the Lagrangian equations for 
x{q,T) near t — 0, which can be linearized because the 
displacement x — q is small, it is easily shown that slaving 
implies the absence of the 'decaying mode,' which behaves as 
T~^l'^ in an Einstein-de Sitter universe and is thus singular 
at r = (for details see Appendix IXjl. 

Slaving is also a sufficient condition for the initial prob- 
lem to be well posed. It is indeed easily shown recursively 
that Q-J^ll admit a solution in the form of a formal Tay- 
lor series in r (a related expansion in volving only potentials 
may be found in lCatelan et alJll995^ : 

v(x,t) = w''^-' (a;) + rw'"""' (a;) + r'^w'^' (a;) + • • • , (5) 
^g(x,r) = 4°)(x)+r4^)(x)+rVf'(a^) + ---, (6) 
p(x,r) = 1 + V(a;)+T^p<^)(x) + ---. (7) 

Furthermore, is easily shown to be curl-free for 

any n. 

Several important consequences of slaving extend to 
later times as 'fossils' of the earliest dynamics. First, as al- 
ready stressed in the Introduction, the whole dynamics is 
determined by only one scalar field (e.g., the initial gravita- 
tional potential) which we can hope to determine from the 
knowledge of the present density field. 

Second, slaving trivially rules out multi-streaming up 
to the time of formation of caustics. Since we are working 
with coUisionless matter, the dynamics should in principle 
be governed by the Vlassov-Poisson'^ kinetic equation which 
allows at each (a;, r) point a non-trivial distribution function 
f{x,v,T). Slaving selects a particular class of solutions for 
which the distribution function is concentrated on a single- 
speed manifold, thereby justifying the use of the Euler- 
Poisson equation without havi ng to invo ke any hydrody- 
namical limit (see, e.g., Vcraasso la et alll9 94: Catcl an et alJ 
[l993). 

Third, it is easily checked from Q that the initial slaved 
velocity, which is obviously curl-free, remains so for all later 
times (up to formation of caustics). Note that this vanishing 
of the curl holds in Eulerian coordinates. A similar property 
in Lagrangian coordinates can only hold approximately but 
will play an important role in the sequel (Section 1^. 



terms of a single-speed solution to the Euler-Poisson equa- 
tions. In A''-body simulations, multi-stream regions are usu- 
ally found to be of relatively small extension in one or several 
space directions, typically not more than a few Mpc, and 
hence have a small volume, although th ey contain a signifi- 
cant fraction of the total mass (see, e.g. Iweinbere fc Guiml 
1 199m . 

In order not to have to deal with tiny multi-stream re- 
gions, we replace the true mass distribution by a 'macro- 
scopic' one which has a regular part and a singular (col- 
lapsed) part, the latter concentrated on objects of dimension 
less than three, such as points or lines. 

The general problem of reconstruction is to find as much 
information as possible on the history of the evolution that 
carries the initial uniform density into the present macro- 
scopic mass distribution, including the evolution of the ve- 
locities. In principle we would like to find a solution of the 
Euler-Poisson initial-value problem leading to the present 
density field po{x). 

A more restricted problem, which we call the 'dis- 
placement reconstruction,' is to find the Lagrangian map 
q i—> x{q) and its inverse x i— » q{x), or in other words to 
answer the question: where does a given 'Monge molecule'^ 
of matter originate from? Of course, the inverse Lagrangian 
map will not be single-valued on mass concentrations. Fur- 
thermore, for practical cosmological applications, we define 
a 'full reconstruction problem' as (i) displacement recon- 
struction and (ii) obtaining the initial and present peculiar 
velocity fields, Vin{q) and vo{x). 

We shall show in this paper that the displacement re- 
construction problem is uniquely solvable and that the full 
reconstruction problem has a unique solution outside of mass 
concentrations; as to the latter, they are traced back to col- 
lapsed regions in the Lagrangian space whose shape and po- 
sitions are well defined but the inner structure of density 
and velocity fiuctuations is irretrievably lost. 



3 POTENTIAL LAGRANGIAN MAPS: THE 
MAK RECONSTRUCTION 

In this and the next two sections we shall assume that the 
Lagrangian map from initial positions to present ones is po- 
tential 

X = V,<E-(r7), (8) 

and furthermore that the potential $(q) is convex, which is, 
as we shall see, related to the absence of multi-streaming. 



2.2 Formulation of the reconstruction problem 

The present Universe is replete with high-density structures: 
clusters (point-like objects), filaments (line-like objects) and 
perhaps sheets or walls.* 

The internal structure of such mass concentrations cer- 
tainly displays multi-streaming and cannot be described in 

^ Actually written for the first time by Ijeanj il91Sl) . 

* Whetiier the Great Wall and the Sculptor Wall are shee t-like 

or filamont-like is a moot point iSathvaprakash et al.lll'99S^ . 



3.1 Approximations leading to maps with convex 
potentials 

The motivation for the pote ntial assumption, first used by 
iBertschinger fc pekel| l| l989^■^ comes from the Zel'dovich 
approximation l|Zerdov'ichlll97ff l. denoted here by ZA, and 
its refinements. To recall how the ZA comes about, let us 

^ For Monge and his contemporaries, the word 'molecule' meant 
a Leibniz infinitesimal element of mass; see Appendix El 
^ In connection with what was called later the La grangian PO- 
TENT method jPekel. Bertschinger fc Fabed[T99cl) . 
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start from the equations for the Lagrangian map x{q,T), 
written in the Lagrangian coordinate q f Appendix IXl 

Dlx = -^{Drx + V^^s), (9) 
Zt 

Vl^s = i [(detV,a;)-i-l] , (10) 

where Dt is the Lagrangian time derivative and = 
{dqj/dxi)Vqj is the Eulerian gradient rewritten in La- 
grangian coordinates. As shown in Appendix^ in one space 
dimension the Hubble drag term DrX and the gravitational 
acceleration term Vcc<^g cancel exactly. Slaving, discussed 
in Section 12.11 means that the same cancellation holds to 
leading order in any dimension for small r. The ZA extends 
this as an approximation without the restriction of small r. 
Within the ZA, the acceleration D^x vanishes. Hence the 
Lagrangian map has the form 

x{q,T) = q + T{Drx)in{q) = q-T\7gip'^\q) (11) 
= V,$((7,r) 

with the potential 

<l.(q,r) = ^-r^-''(q). (12) 

Furthermore, taking the time derivative of (IIH . we see that 
the velocity DTx{q,T) is curl-free with respect to the La- 
grangian coordinate q. 

Potentiality of the Lagrangian map (and consequently 
the Lagrangian potentiality of the velocity) is perhaps the 
most important feature of the ZA. Unlike the vanishing of 
the acceleration, it does not depend on the choice of the 
linear growth factor as the time variable. However, unac- 
celerated but vortical flow would fail to exhibit the can- 
cellation necessary for the ZA to hold. It is noteworthy 
that the potentiality is not limited to the ZA: indeed, the 
latter can be formulated as the first order of a system- 
atic Lagrangian perturbation theory in which, up to sec- 
ond order, the Lagra ngian map is still potential under slav- 
ing (Moutar de et aljrifl 91; Buchcrt 1992; Buchcrt fc Ehlerj 
ll99,Sli : .Munshi. Sahni fc Starobinskvj.1994 : .Catelan.l995i' l . 

It is well known that the ZA map defined by IIIH ceases 
in general to be invertible due to the formation of multi- 
stream regions bounded by caustics. Since particles move 
along straight lines in the ZA, the formation of caustics 
proceeds just as in ordinary optics in a uniform medium 
in which light rays are also straight.^ One of the prob- 
lems with the ZA is that caustics, which start as localized 
objects, quickly grow in size and give unrealistically large 
multi-stream regions. 

A modification of the ZA that has no multi-streaming 
at all, but sharp mass concentrations in the form of 
shocks and other s i ngular ities, has been introduced by 
Gurbatov fc Saichev| l|l 984|: s ee also iGurbatov et alJ Il989l : 
Shandarin fc Ze^dovichlll989^ . It is known as the adhesion 
model. In Eulerian coordinates it amounts to using a multi- 
dimensional Burgers equation (see, e.g..|Frisch & Bee 20o3) 

OtV + {V ■ \7a:)v = i^V^W, V = —V^fv, (13) 

^ Catastrophe theory has been used to classify the different types 
of sin gularities thus obtained lArnoI'd. Shandarin fc Zel'dovichl 
Il983). 



taken in the limit where the viscosity v tends to zero. In 
Lagrangian coordinates, the adhesion model is obtained 
from the ZA by replacing the velocity potential ^{q,t) 
given by 11211 by its co nvex hull $c{q,t) in the q variable 
JVergassola et alJll99J) . 

Convexity is a concept which plays an important role 
in this paper, and a few words on it are in order here (see 
also Appendix ICH . A body in the three-dimensional space 
is said to be convex if, whenever it contains two points, it 
contains also the whole segment joining them. A function 
f{q) is said to be convex if the set of all points lying above 
its graph is convex. The convex hull of the function $(q) 
is defined eis the largest convex function whose graph lies 
below that of $(q). In two dimensions it can be visualized 
by wrapping the graph of $(9) tightly from below with an 
elastic sheet. 

Note that "l>(q,r) given by 1121 is obviously convex for 
small enough r since it is then very close to the parabolic 
function |<7p/2. After caustics form, convexity is lost in the 
ZA but recovered with the adhesion model. It may then 
be shown that those regions in the Lagrangian space where 
$(q, t) does not coincide with its convex hull will be mapped 
in the Eulerian space to sheets, lines and points, each of 
which contains a finite amount of mass. At these locations 
the Lagrangian map does not have a uniquely defined La- 
grangian antecedent but such points form a set of vanishing 
volume. Everywhere else, there is a unique antecedent and 
hence no multi-streaming. 

Although the adhesion model has a number of 
known shortcomings, such as non-conservation of momen- 
tum in more than one dimension, it has been found 
to be in better agreemen t with A^-body simulations 
than the ZA (j^einbcre fc GunnI [i§90). Other single- 
speed approximations to multi-stream fiow, overcoming 
difficultie s of the adhesion model, are discussed e.g. by 
Shandarin_j^_ggjjr^gj^rg^asb|j[l 9 : iBuchert fc Domineue j 
(11993); IFanelli fc Aurelll J2002ir in such models, multi- 
streaming is completely suppressed by a mechanism of mo- 
mentum exchange between neighbouring streams with differ- 
ent velocities. This is of course a common phenomenon in or- 
dinary fluids, where it is due to viscous diffusion; dark mat- 
ter is however essentially coUisionless and the usual mech- 
anism for generating viscosity does not operate, so that 
a non-coUisional mechanism must be invoked. A qualita- 
tive explanation using the modification of the gravitational 
forces after the formation of ca ustics has been proposed by 
IShandarin fc Zel'dovichl if 1989'). In our opinion the mecha- 
nism limiting multi-streaming to rather narrow regions is 
poorly understood and deserves considerable further inves- 
tigation. 

3.2 The Monge— Ampere equation: a consequence 
of mass conservation and potentiality 

We now show that the assumption that the Lagrangian map 
is derived from a convex potential leads to a pair of non- 
linear partial differential equations, one for this potential 
and another for its Legendre transform. 

Let us flrst assume that the present distribution of 
mass has no singular part, an assumption which we shall 
relax later. Since in our notation the initial quasi-uniform 
mass distribution has unit density, mass conservation im- 
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plies po(a;)d'^a; = d'^q, which can be rewritten in terms of 



detVqa; = — , ^ . (14) 

Under the potential assumption (0, this takes the form 
1 



det(V,.V„.$((7)) = 



po(V,$(q)) 



(15) 



A similar equation follows also from Eqs. (1) and (2) of 
iBertschinger fc DekeJ (|l98fll ). 

A simpler equation, in which the unknown appears only 
in the left-hand side, viz Eq. (1191 below, is obtained for the 
potential of the inverse Lagrangian map q{x). Key is the 
observation that the inverse of a map with a convex poten- 
tial has also a convex potential, and that the two potentials 
are Legendre transforms of each other.* A purely local proof 
of this statement is to observe that potentiality of q{x) is 
equivalent to the symmetry of the inverse Jacobian matrix 
Vxq which follows because it is the inverse of the symmet- 
rical matrix VqX; convexity is equivalent to the positive- 
definiteness of these matrices. Obviously the function 



6(a;) = X ■ q{x) ~ <l>(q(a;)). 



(16) 



which is the Legendre transform of ^{q), is the potential 
for the inverse Lagrangian map. The modern definition of 
the Legendre transformation (see Appendix IClfl . needed for 
generalization to non-smooth mass distributions, is 



0(a;) = max a; ■ q — ^(q), 
<l?(q) = max x ■ q — 0{x). 



(17) 
(18) 



In terms of the potential O, mass conservation is imme- 
diately written as 



det(V,,V,^e(a;)) = po{x). 



(19) 



This equation, which has the determinant of the second 
derivatives of the unknown in the left-hand side and a pre- 
scribed (positive) function in the right-hand side, is called 
the (elliptic) Monge-Ampere equation (see Appendix IbI for 
a historical perspective). 

Notice that our Monge-Ampere equation may be 
viewed as a non-linear generali zation of the Po i sson e qua- 
tion (used for reconstruction bv lNusser fc Dekell l)l992h : see 
also Section lY.lH . to which it reduces if particles have moved 
very little from their initial positions. 

In actual reconstructions we have to deal with mass con- 
centration in the present distribution of matter. Thus the 
density in the right-hand side of 11911 has a singular com- 
ponent (a Dirac distribution concentrated on sets carrying 
the concentrated mass) and the potential Q ceases to be 
smooth. As we now show, a generalized meaning can never- 
theless be given to the Monge-Ampere equation by using the 



° Besides our problem, this fact prominently appears in two 
other fields of physics: in classical mechanics, the Lagrangian and 
Hamiltonian functions are Legendre transforms of each other - 
their gradients relate the generalized velocity and momentum - 
and so are, in thermodynamics, the internal energy and the Gibbs 
potential, implying the same relation between extensive and in- 
tensive parameters of state. 



key ingredient in its derivation, namely mass conservation, 
in integrated form. 

For a nonsmooth convex potential O, taking the gra- 
dient Va,0(x) still makes sense if one allows it to be mul- 
tivalued at points where the potential is not differentiable. 
The gradient at such a point x is then the set of all pos- 
sible slopes of planes touching the graph of Q at {x,Q{x)) 
(this idea is given a precise mathematical formulation in Ap- 
pendix IcT^ . As x varies over an arbitrary domain De in the 
Eulerian space, its image q{x) sweeps a domain q{T>^) in 
the Lagrangian space, and mass conservation requires that 



po{x) d a; = 



A'q, 



(20) 



where we take into account that q{x) — Va;0(x). Eq. 1201 
must hold for any Eulerian domain T>e; this requirement 
is known as the weak formulation of the Monge-Ampere 
equation 11911 . A symmetric formulation may be written for 
list in terms of x{q) = Vo 'I'(q). For further material on the 
weak formulation see, e.g.. IPoeorelovl (|l978h . 

Considerable literature has been devote d to the Monge- 
Amporo equation in recent years (see, e.g.. lCaffareiiilll99gl : 
Caffarclli & Milman 1999). We mention now a few results 
which are of direct relevance for the reconstruction problem. 

In a nutshell, one can prove that when the domains oc- 
cupied by the mass initially and at present are bounded and 
convex, the Monge-Ampere equation - in its weak formu- 
lation - is guaranteed to have a unique solution, which is 
smooth unless one or both of the mass distributions is non- 
smooth. The actual construction of this solution can be done 
by a variational method discussed in the next section. 

A similar result holds also when the present density field 
is periodic and the same periodicity is assumed for the map. 

Als o relevant, as we shall see in Section [3.41 is a recent 
result of lCaffarelli fc Lil (|2Q^) : if the Monge-Ampere equa- 
tion is considered in the whole space, but the present density 
contrast S — p — 1 vanishes outside of a bounded set, then 
the solution 0(a;) is determined uniquely up to prescription 
of its asymptotic behaviour at infinity, which is specified by 
a quadratic function of the form 



(x) = {x, Ax) + {b, x) + c. 



(21) 



for some positive definite symmetric matrix A with unit de- 
terminant, vector b and constant c. 



3.3 Optimal mass transportation 

As we are going to see now, the Monge-Ampere equation 
l|19|l is equivalent to an instance of what is called the 'opti- 
mal mass transportation problem.' Suppose we are given two 
distributions pin(q) and po(x) of the same amount of mass 
in two three-dimensional convex bounded domains T>i^ and 
Do. The optimal mass transportation problem is then to find 
the most cost-effective way of rearranging by a suitable map 
one distribution into the other, the cost of transporting a 
unit of mass from a position q G Din to a; £ Do being a 
prescribed function c(q,x). 

Denoting the map by x{q) and its inverse q{x), we can 
write the problem as the requirement that the cost 

1= c{q,x{q))pi^{q)A\= / c{q{x) , x) po{x) A'' X (22) 



© 0000 RAS, MNRAS 000, 000-000 



6 Y. Brenier et al. 



be minimum, with the constraints of prescribed 'terminal' 
densities pin and po and of mass conservation pin(<j')d"^q = 
po{x)A'x.'> 

This problem goes back to iMongd il78 J) who consid- 
ered the case of a linear cost function c(q, x) = \x — q\ (see 
Appendix IbI and Fig.^. 

For our purposes, the central result is that the problem 
of finding a potential Lagrangian map with presecribed initial 
and present mass density fields is equivalent to a mass trans- 
ortation problem w ith quadratic cost. Indeed, it is known 
Brenieill987l . ll99ll) that, when the cost is a quadratic func- 
tion of the distance, so that 



!^^^po(a;)d^a;,(23) 



the solution q{x) to the optimal mass transportation prob- 
lem is the gradient of a convex function, which then must 
satisfy the Monge-Ampere equation H19|l by mass conserva- 
tion. 

A particularly simple variational proof can be given for 
the smooth case, when the two mutually inverse maps x{q) 
and q{x) are both well defined. 

Performing a variation of the map £c(q), we cause a 
mass element in the Eulerian space that was located at 
x{q) to move to x{q) + Sx{q). This variation is constrained 
not to change the density field po. To express this con- 
straint it is convenient to rewrite the displacement in Eule- 
rian coordinate Sxe{x) = Sx{q{x)). Noting that the point 
X gets displaced into y = x + Sx, we thus require that 
po{x) d^x = po{y) d^y or 

poix) = po{x + 5xe{x)) det(Va;(a; + 5a;E(a;))). (24) 
Expanding this equation, we find that, to the leading order. 



{po{x) Sxe{x)) = 0, 



(25) 



an equation which just expresses the physically obvious fact 
that the mass flux po{x) Sxe{x) should have zero divergence. 
Performing the variation on the functional / given by (I23II . 
we get 



SI 



{x{q) - q) ■ Sx{q) pin{q) d q 



/ {x - q{x)) ■ (^po{x) Sxe{x)) d^x = 0, 

J-Dn 



(26) 



which has to hold under the constraint 1251 1. In other words, 
the displacement x — q{x) has to be orthogonal (in the 
L2 functional sense) to all divergence-less vector fields and, 
thus, must be a gradient. Since x is obviously a gradient, it 
follows that q{x) = \/xQ{x) for a suitable potential Q. 

It remains to prove the convexity of Q. First we prove 
that the map x i—> q{x) = Va;0(a;) is monotone, i.e., by 
definition, that for any xi and X2 



(x2 - xi) ■ {q{x2) - q{xi)) ^ 0. 



(27) 



Indeed, should this inequality be violated for some xi,X2, 
the continuity of q{x) would imply that for all xi,X2 close 
enough to Xi , X2 



\q{xi) - Xi\'^ + \q(x2) - X2\^ 

> \q{x2) - Xi\"' + \q{xi) - X2f ■ 



(28) 



This in turn means that if we interchange the destinations 
of small patches around xi and X2, sending them not to 
the corresponding patches around q{x\) and q(x2) but vice 
versa, then the value of the functional I will decrease by 
a small yet positive quantity, and therefore it cannot be 
minimum for the original map.^'^ 

To complete the argument, observe that convexity of 
a smooth function Q{x) follows if the matrix of its second 
derivatives V^i^ xj Q{x) is positive definite for all x. Substi- 
tuting q{x) — \7a:Q{x) into 1271 1. assuming that X2 is close 
to xi and Taylor expanding, we find that 



{x2 - xi) ■ (Vx.Va;^6(a;i) (x2 - xi)) ^ 0. 



(29) 



As X2 is arbitrary, this proves the desired positive definite- 
ness and thus establishes the equivalence of the Monge- 
Ampere equation H19|l and of the mass transportation prob- 
lem with quadratic cost. 

This equivalence is actually proved unde r much weaker 
conditions, not requiring any smoothness ("Brenier 'l987|, 
[1991). The proof makes use of the 'relaxed' rcfornmlation of 
the mass transportation problem due to lKantorovichI lll943) . 
Instead of solving the highly non-linear problem of finding 
a map q{x) minimizing the cost H22|l with prescribed ter- 
minal densities, Kantorovich considered the linear program- 
ming problem of minimizing 



1=1 / c{q,x)p{q,x)d'^qd^x, 



(30) 



under the constraint that the joint distribution p(q, a;) is 
nonnegative and has marginals pin(<3') and po{x), the latter 
being equivalent to 



p{q,x)d X = pi^{q), I p{q,x) d q = po{x). (31) 

■Co •''DiTi 

Note that if we assume any of the two following forms for 
the joint distribution 

p{q,x) = po{x)5{q - q{x)) 
p{q, x) = pin(q) 5{x- x{q)) , 



(32) 



we find that / reduces to the cost I as defined in (j22j. This 
relaxed formulation allowed Kantorovich to establish the ex- 
istence of a mimimizing joint distribution. 

The relaxed formulation can be used to show that the 
minimizing solution actually defines a map, which need not 
be smooth if one or both of the terminal distribution have a 
singular component (in our c ase, when mass con centrations 
are present). The derivation iBrenieijll987l . Il99 J) makes use 
of the technique of duality (Appendix IC2II . which will also 
appear in discussing algorithms (Section I4.2|l and recon- 
struction beyond the potential hypothesis (Section (SJ. 

We have thus shown that the Monge-Kantorovich opti- 
mal mass transportation problem can be applied to solving 
the Monge-Ampere equation. The actual implementation 



^ Note that x(q) = q does not solve the above problem as it vio- As we shall see in Section |4. II the converse is not true: mono- 

lates the latter constraint unless the terminal densities are iden- tonicity alone does not imply that the integral / is a minimum; 

tical. the minimizing map must also be potential. 
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Figure 2. A one-dimensional example of non-unique reconstruc- 
tion of the Lagrangian map in the presence of multi-streaming. 
The density distribution (upper graph) is generated by a multi- 
streaming Lagrangian map (thick line of lower graph) but may 
also be generated by a spurious single-stream Lagrangian map 
(dashed line). 

(Sectional, done for a suitable discretization, will be hence- 
forth called Monge-Ampere-Kantorovich (MAK). 

3.4 Sources of uncertainty in reconstruction 

In this section we discuss various sources of non-uniqueness 
of the MAK reconstruction: multi-streaming, collapsed re- 
gions, reconstruction from a finite patch of the Universe. 

We have stated before that our uniqueness result applies 
only in so far as we can treat present-epoch high-density 
multi-stream regions as if they were truly collapsed, ignor- 
ing their width. We now give a simple one-dimensional ex- 
ample of non-uniqueness in which a thick region of multi- 
streaming is present. Fig.|5|shows a multi-stream Lagrangian 
map and the associated density distribution; the inverse 
map qix) is clearly multi-valued. The same density distribu- 
tion may however be generated by a spurious single-stream 
Lagrangian map shown on the same figure. There is no way 
to distinguish between the two inverse Lagrangian maps if 
the various streams cannot be disentangled. 

Suppose now that the present density has a singular 
part, i.e. there are mass concentrations present which have 
vanishing (Eulerian) volumes but possess finite masses. Ob- 
viously any such object originates from a domain in the 
Lagrangian space which occupies a finite volume. A one- 
dimensional example is again helpful. Fig. shows a La- 
grangian map in which a whole Lagrangian shock interval 
[51,52] has collapsed into a single point of the x axis. Out- 
side of this point the Lagrangian map is uniquely invertible 
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Figure 3. Two initial velocity profiles v\^{q) (bottom, solid and 
dashed lines) leading to the same Lagrangian map x = q + TV\^{q) 
(top, solid line) in the adhesion approximation. The Zel'dovich 
approximation would give multistreaming (top, dashed line). 
Hatched areas (bottom) are equal in the adhesion dynamics. 

but the point itself has many antecedents. Note that the 
graph of the Lagrangian map may be inverted by just inter- 
changing the q and x axes, but its inverse contains a piece 
of vertical line. The position of the Lagrangian shock in- 
terval which has collapsed by the present epoch is uniquely 
defined by the present mass field but the initial velocity fiuc- 
tuations in this interval cannot be uniquely reconstructed. 
In particular there is no way to know if collapse has started 
before the present epoch. We can of course arbitrarily as- 
sume that collapse has just happened at the present epoch; 
if we also suppose that particles have travelled with a con- 
stant speed, i.e. use the Zel'dovich/adhesion approximation, 
then the initial velocity profile within the Lagrangian shock 
interval will be linear (Fig. |^ . Any other smooth velocity 
profile joining the same end points would have points where 
its slope (velocity gradient) is more negative than that of the 
linear profile (Fig. (SJ and thus would have started collapse 
before the present epoch (in one dimension caustics appear 
at the time which is minus the inverse of the most negative 
initial velocity gradient). 

All this carries over to more than one dimension. The 
MAK reconstruction gives a unique antecedent for any 
Eulerian position outside mass concentrations. Each mass 
concentration in the Eulerian space, taken globally, has a 
uniquely defined Lagrangian antecedent region but the ini- 
tial velocity field inside the latter is unknown. In other 
words, displacement reconstruction is well defined but full 
reconstruction, based on the Zel'dovich/adhesion approxi- 
mation for velocities, is possible only outside of mass concen- 
trations (note however that velocities in the Eulerian space 



© 0000 RAS, MNRAS 000, 000-000 



8 Y. Brenier et al. 



are still reconstructed at almost all points). We call the cor- 
responding initial Lagrangian domains collapsed regions. 

Finally, we consider a uniqueness problem arising from 
knowing the present mass distribution only truncated over 
a finite Eulerian domain Do , as is necessarily the case when 
working with a real catalogue. If we also know the cor- 
responding Lagrangian domain I?in and both domains are 
bounded and convex, then uniqueness is guaranteed (see Sec- 
tion EiU. What we know for sure about Din is its volume, 
which (in our units) is equal to the total mass contained 
in Do- Its shape and position may however be constrained 
by further information. For example, if we know that the 
typical displacement of mass elements since decoupling is 
about ten Mpc in comoving coordinates (see Section]^ and 
our data extend over a patch of typical size one hundred 
Mpc, then there is not more than a ten percent uncertainty 
on the shape of Din. Additional information about peculiar 
velocities may also be used to constrain Din. 

Note also that a finite-size patch Do with unknown an- 
tecedent Din will give rise to a unique reconstruction (up to 
a translation) if we assume that it is surrounded by a uni- 
form backgrou nd extending to infinity. This is a consequence 
of the result of lCaffarelli fc lI mentioned at the end of Sec- 
tion l.'i.2l The arbitrary linear term in l|21|l corresponds to 
a translation; as to the quadratic term, it is constrained by 
the cosmological principle of isotropy to be exactly |q|^/2. 



4 THE MAK METHOD: DISCRETIZATION 
AND ALGORITHMICS 

In this section we show how to compute the solution to the 
Monge-Ampere-Kantorovich (MAK) problem the known 
present density field. First the problem is discretized into 
an assignment problem (Section 14.1^ . then we present some 
general tools which make the assignment problem computa- 
tionally tractable ('Section l4.2ll and finally we present, to the 
best of our knowledge, the most effective method for solving 
our particular assignment problem, based on the auction 
algorithm of D. Bertsekas (Section 14. 3L and details of its 
implementation for the MAK reconstruction (Section 14.4^ . 

4.1 Reduction to an assignment problem 

Perhaps the most natural way of discretizing a spatial mass 
distribution is to approximate it by a finite system of identi- 
cal Dirac point masses, with possibly more than one mass at 
a given location. This is compatible both with A'^-body simu- 
lations and with the intrinsically discrete nature of observed 
luminous matter. Assuming that we have A'^ unit masses 
both in the Lagrangian and the Eulerian space, we may write 

N N 

pa{x) ^'^5{x - Xi), pin{q) ^^5{q - Qj). (33) 

For discrete densities of this form, the mass conservation 
constraint in the optimal mass transportation problem (Sec- 
tion I3.3|l requires that the map q{x) induce a one-to-one 
pairing between positions of the unit masses in the x and q 
spaces, which may be written as a permutation of indices 
that sends Xi t o g, ;,;). Substituting this into the quadratic 
cost functional 1231 . we get 



iV I |2 

1=1 

We thus reduced the problem to the purely combinatorial 
one of finding a permutation (or its inverse i{j)) that 
minimizes the quadratic cost function 1341 . 

This problem is an instance of the general assignment 
problem in combinatorial optimization: for a cost matrix dj, 
find a permutation j{i) that minimizes the cost function 

JV 

7 = ^Cij(i). (35) 
1=1 

As we shall see in the next sections, there exist effective 
algorithms for finding minimizing permutations. 

Before proceeding with the assignment problem, we 
should mention an alternative approach in which discretiza- 
tion is performed only in the Eulerian space and the ini- 
tial mass distribution is kept continuous and uniform. Min- 
imization of the quadratic cost function will then give rise 
to a tesselation of the Lagrangian space into polyhedric re- 
gions which end up collapsed into the discrete Eulerian Dirac 
masses. Basically, the reason why these regions are poly- 
hedra is that the convex potential ^{q) of the Lagrangian 
map has a gradient which takes only finitely many values. 
This problem, which has been studie d by Aleksandrov and 
Pogorelov (see, e.g., IPoeorelovl Il978f) . is closely related to 
Minkowski's (1897.) famous problem of constructing a con- 
vex polyhedron with prescribed areas and orientations of 
its faces (in our setting, areas and orientations correspond 
to masses and values of the gradient). Uniqueness in the 
Minkowski problem is guaranteed up to a translation. Start- 
ing with Minkowski's own very elegant solution, various 
methods of constructing solutions to such geometrical ques- 
tions have been devised. So far, we have not been able to 
make use of such ideas in a way truly competitive with dis- 
cretization in both spaces and solving then the assignment 
problem. 

The solution to our assignment problem (with quadratic 
cost) has the important property that it is monotone: for 
any two Lagrangian positions q^ and the corresponding 
Eulerian positions xi and X2 are such that 

(cci - 0^2) ■ (Qi - Qa) > 0. (36) 

This is of course the discrete counterpart of (1271 . In one di- 
mension, when all the Dirac masses are on the same line, 
monotonicity implies that the leftmost Lagrangian position 
goes to the leftmost Eulerian position, the second leftmost 
Lagrangian position to the second leftmost Eulerian posi- 
tion, etc. It is easily checked that this correspondence mini- 
mizes the cost 1341 . 

In more than one dimension, a correspondence between 
Lagrangian and Eulerian positions that is just monotone 
will usually not minimize the cost (a simple two-dimensional 
counterexample is given in Fig. |1J.^^ Actually, a much 
stronger condition, called cyclic monotonicity, is needed in 
order to minimize the cost. It requires fc- monotonicity for 

Note that in one dimension, in the continuous case, any map 
is a gradient and we have already observed in Section IS . 3l that if a 
gradient map is monotone it is the gradient of a convex function. 
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Figure 4. Two monotone assignments sending white points to 
black ones: (a) an assignment that is vastly non-optimal in terms 
of quadratic cost but cannot be improved by any pair interchange; 
(b) the optimal assignment, shown for comparison. 



any k between 2 and A'^; the latter is defined by taking any 
k Eulerian positions with their corresponding Lagrangian 
antecedents and requiring that the cost 1341 1 should not de- 
crease under an arbitrary reassignment of the Lagrangian 
positions within the set of Eulerian positions taken. Note 
that the usual monotonicity corresponds to 2-monotonicity 
(stability with respect to pair exchanges). 

A strategy called PIZA (Path Interchange Zel'dovich 
Approximation) for constructing monotone correspondences 
between Lagrangian and Eulerian positions has been pro- 
posed bv lCroft fc Gaztafiagal Il997f) . In PIZA, a randomly 
chosen tentative correspondence between initial and final 
positions is successively improved by swapping randomly 
selected pairs of initial particles whenever H36^ is not sat- 
isfied. After the cost 1341 ceases to decrease between iter- 
ations, an approximation to a monotone correspondence is 
established, which is generally neither unique, as already 
observed bv lValentine. Saunders fc Tavloil i200Ci) in testing 
PIZA reconstruction, nor optimal. We shall come back to 
this in Sections 13 and 



4.2 Nuts and bolts of solving the assignment 
problem 

For a general set of A'^ unit masses, the assignment prob- 
lem with the cost function H34^ has a single solution which 
can obviously be found by examining all A^! permutations. 
However, unlike computationally hard problems, such as the 
travelling salesman's, the assignment problem can be han- 
dled in 'polynomial time' - actually in not more than 0{N'^) 
operations. All methods achieving this use a so-called dual 
formulation of the problem, based on a relaxation similar to 
that applied by Kantorovich to the optimal mass transporta- 
tion fSection 13.31 a brief introduction to duality is given in 
Appendix IC2^ . In this section we explain the basics of this 
technique, using a variant of a simple mechanical model in- 
troduced in a more general setting bv lHenonI (^995, 2002). 

Consider the general assignment problem of minimizing 
the cost 1351 1 over all permutations We replace it by a 
'relaxed,' linear programming problem of minimizing 



(37) 




columns 



Figure 5. An analogue computer solving the assignment problem 
for Af = 4. 



(38) 



for all i, j, an obvious discrete analogue of H31^ . We show 
now that it is possible to build a simple mechanical device 
(Fig.|KJ which solves this relaxed problem and that the solu- 
tion will in fact determine a minimizing permutation in the 
original assignment problem (i.e., for any i or j fixed, only 
one fij will be unit and all other zero). The device acts as 
an analogue computer: the numbers involved in the problem 
are represented by physical quantities, and the equations are 
replaced by physical laws. 

Define coordinate axes x, y, z in space, with the z axis 
vertical. We take two systems of A'' horizontal rods, parallel 
to the X and y axes respectively, and call them columns and 
rows, referring to columns and rows of the cost matrix. Each 
rod is constrained to move in a corresponding vertical plane 
while preserving the horizontal orientation in space. For a 
row rod Ai, we denote the z coordinate of its bottom face 
by Qi and for a column rod Bj, we denote the z coordinate 
of its top face Pj. Row rods are placed above column rods, 
therefore Ui ^ /3j for all i,j (see Fig.jKJ. 

Upper (row) rods are assumed to have unit weight, 
and lower (column) rods to have negative unit weight, or 
unit 'buoyancy' Therefore both groups of rods are subject 
to gravitational forces pulling them together. However, this 
movement is obstructed by A'^^ small vertical studs of negli- 
gible weight put on column rods just below row rods. A stud 
placed at projected intersection of column Bj and row Ai 



has length C - 



with a suitably large positive constant C 



where auxiliary variables fij satisfy 



and thus constrains the quantities ai and f3j to satisfy the 
stronger inequality 

a,-l3j^C-c,j. (39) 
The potential energy of the system is, up to a constant. 
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(40) 



3=1 



In linear programming, the problem of minimizing 1401 un- 
der the set of constraints given by 13911 is called the dual 
problem to the 'relaxed' one 

(EB-lEHJ (see Appendix OJ; 
the a and (3 variables are called the dual variables. 

The analogue computer does in fact solve the dual prob- 
lem. Indeed, first hold the two groups of rods separated from 
each other and then release them, so that the system starts 
to evolve. Rows will go down, columns will come up, and 
contacts will be made with the studs. Aggregates of rows 
and columns will be progressively formed and modified as 
new contacts are made, giving rise to a complex evolution. 
Eventually the system reaches an equilibrium, in which its 
potential en ergv l|40t is m inimum and all constraints 13911 
are satisfied (|Henonll2002t l. Moreover, it may be shown that 
the solution to the original problem 137|I - H38^ is expressible 
in terms of the forces exerted by the rods on each other 
at equilibrium and is typically a one-to-one correspondence 
between the AiS and the BjS (for details, see ADDendix lC3|l . 

The common feature of many existing algorithms for 
solving the assignment problem, which makes them more 
effective computationally than the simple enumeration of 
all A^! permutations, is the use of the intrinsically contin- 
uous, geometric formulation in terms of the pair of linear 
programming problems |(23~(EHJ a-nd I)40|l - lj39|l . The me- 
chanical device provides a concrete model for this formula- 
tion; in fact, assignment algorithms can be regarded as de- 
scriptions of specific procedures to make the machine reach 
its equilibrium state. An introduction into algorithmic as- 
pects of solving the assignment problem, including a proof of 
the 0{N^) theoretical bound on the number of operations, 
based o n the Hungarian method of .Kuhn ( 1955i) . may be 
found in IPaoadimitrioxi fc Steiglit3 l)l982|) 7 

In spite of the general 0{N^) theoretical bound, various 
algorithms may show very different performance when ap- 
plied to a specific optimization problem. Duringthe prepa- 
ration of the earli er publication I Frisch et al...,2002f) the dual 
simplex method of lBalins'kil il986t) wa s used, with so me mod- 
ifications inspired by algorithm B of iHenonI (|22q2)- Several 
other algorithms were tried subsequently, including an adap- 
tation o f algorithm A of the latt er reference and the algo- 



rithm of 
work of 



^jrk^r^ 



Tbmizaw 



Dgrig^ ()l98ff l , itself based on the earlier 
lll97lh . For the time being, the fastest 
running code by far is based on the auction algorithm of 
[Scrtsckas (1992, 2001), arguably the most effective of ex- 
isting ones, which is discussed in the next section. Needless 
to say, all these algorithms arrive at the same solution to 
the assignment problem with given data but can differ by 
several orders of magnitude in the time it takes to complete 
the computation. 



4.3 The auction algorithm 

We explain here the essense of the auction algorithm in 
terms of our mechanical devi ce. Note that the original p re- 
sentation of this algorithm ljBertsekajll98ll Il99l 120011^ is 
based on a different perspective, that of an auction, in which 
the optimal assignment appears as an economic rather than 
a mechanical equilibrium; the interested reader will benefit 
much from reading these papers. 

Put initially the column rods at zero height and all row 
rods well above them, so that no contacts are made and con- 
straints H39|l are satisfied. To decrease the potential energy, 
let now the row rods descend while keeping the column rods 
fixed. Eventually all row rods will meet studs placed on col- 
umn rods and stop. Some column rods may then come in 
contact with multiple row rods. Such rods are overloaded: if 
they were not prevented from moving they would descend. 

Note that at this stage any column rod Ai has es- 
tablished a contact with a row rod Bj for which the stud 
length C — Cij is the maximum and the cost Cij the mini- 
mum among other Bs; for Cij = \xi — qj\'^/2, this means that 
any Eulerian position Xi is coupled to its nearest Lagrangian 
neighbour . This coupling is a reasonable guess for the op- 
timal assignment; should it happen to be one-to-one, then 
the equilibrium, and with it the optimal assignment, would 
be reached. It is usually not, so there are overloaded B rods 
and the following procedure is applied to find a compromise 
between minimization of the total cost and the requirement 
of one-to-one correspondence. 

Take any overloaded rod Bj and let it descend while 
keeping other column rods fixed. As Bj descends, row rods 
touching it will follow its motion until they meet studs of 
other column rods and stay behind. The downward motion 
of Bj is stopped only when the last row rod touching Bj is 
about to lose its contact. We then turn to any other over- 
loaded column rod and repeat the procedure as often as 
needed. 

This general step can be viewed as an auction in which 
row rods bid for the descending column rod, offering prices 
equal to decreases in their potential energy as they follow 
its way down. As the column rod descends, thereby increas- 
ing its price, the auction is won by the row rod able to 
offer the largest bidding increment, i.e., to decrease its po- 
tential energy by the largest amount while not violating the 
constraints posed by studs of the rest of column rods. For 
computational purposes it suffices to compute bidding in- 
crements for all competing row rods from the dual a and (3 
variables and assign the descending column rod Bj to the 
highest bidder Ai, decreasing their heights Pj and ai corre- 
spondingly. 

Observe that, at each step, the total potential energy U 
defined by 1401 decreases by the largest amount that can be 
achieved by moving the descending column rod without vio- 
lating the constraints.^* Since l)40|l is obviously nonnegative. 



^■^ This applies to algorithms that never violate constraints 1391 
represented by studs; all practical assignment algorithms known 
to us fall within this category. 



A movie illustrating the subsequent discussion may be found 
at http://www.obs-nice.fr/etc7/movio.html (requires fast Inter- 
net access). 

^"^ This idea of moving a rod, or adjusting a dual variable, up to 
the last point compatible with all the constraints, may be actu- 
ally implemented in a number of ways, giving rise to several pos- 
sible flavours of the auction algorithm. For example, the above 
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the descent cannot proceed indefinitely, and the process may 
be expected to converge quite fast to a one-to-one pairing 
that solves the assignment pr oblem. 

However, as observed bv lBertseka^ l)l98llll99l l20nJ) . 
this 'naive' auction algorithm may end up in an infinite cy- 
cle if several row rods bid for a few equally favourable col- 
umn rods, having thus zero bidding increments. To break 
such cycles and also to accelerate convergence, a perturba- 
tion mechanism is introduced in the algorithm. Namely, the 
constraints l)H9|l are replaced by weaker ones 

a, - /3j > C - c,j - e (41) 

for a small positive quantity e, and in each auction the de- 
scending column rod is pushed down by e in addition to 
decreasing its height by the bidding increment. It can be 
shown that this reformulated process terminates in a finite 
number or rounds; moreover, if all stud lengths are integer 
and e is smaller than 1/A'^, then the algorithm terminates at 
an assignment that is optimal in the unperturbed problem 
(^Bertsekaail99a' l. 

The third ingredient in the Bertsekas algorithm is the 
idea of e-scaling. When the values of dual variables are al- 
ready close to the solution of the dual problem, it usually 
takes relatively few rounds of auction to converge to a solu- 
tion. Thus one can start with large e to compute a rough ap- 
proximation for dual variables fast, without worrying about 
the quality of the assignment, and then proceed reducing e 
in geometric progression until it passes the threshold, 
assuring that the assignment thus achieved solves the initial 
problem. 

Bertsekas' algorithm is especially fast for sparse assign- 
ment problems, in which rods Ai and Bj can be matched 
only if the pair belongs to a given subset A of the set 
of A'^^ possible pairs. We call such pairs valid and define the 
filling factor to be the proportion of valid pairs / = |.4|/A''^. 
When this factor is small, computation can be considerably 
faster: to find the bidding increment for a rod Ai, we need 
only to run over the list of rods Bj such that {i,j) is a valid 
pair. 

Note also that the decentralized structure of the al- 
gorithm facilitates i ts parallelization (see references in 
lBertsekajll993 . l200J) . 

4.4 The auction algorithm for the MAK 
reconstruction 

We now describe the adaptation of the auction algorithm 
to the MAK reconstruction. Experiments with various 
programs contained in Bertsekas' publicly available package 
: / / web . mit . edu /dimitrib / www / auction . txt [ showed 
that the most effective for our problem is AUCTION_flp. It 
assumes integer costs Cij, which in our case requires proper 
scaling of the cost matrix. To achieve this, the unit of length 
is adjusted so that the size of the reconstruction patch 
equals 100, and then the square of the distance between an 

procedure in its most effective implementation requires a parallel 
computer so that groups of several rods can be tracked simultane- 
ously. On sequential computers another, less intuitive procedure, 
in which upper rods are dropped once at a time, proves more 
effective iBertseka£lll992l) . 
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Figure 6. Computing time for different algorithms as a function 
of the numbe r N of points (div i ded b y for normalization). 
Asterisks, the iBurkard fc Derigsl jlflSCl) algorithm (BD); crosses 
and points, the dense and sparse versions of the auction algorithm 
(described in the text). 



initial and a final position is rounded off to an integer. In 
our application, row and column rods correspond to Eule- 
rian and Lagrangian positions, respectively. As the MAK 
reconstruction is planned for application to catalogues of 
10^ and more galaxies, we do not store the cost matrix, 
which would require an 0{N'^) storage space, but rather 
compute its elements on demand from the coordinates, 
which requires only 0{N) space. 

Our problem is naturally adapted for a sparse descrip- 
tion if galaxies travel only a short distance compared to the 
dimensions of the reconstruction patch. For instance, in the 
simulation discussed in Section |3 the r.m.s. distance trav- 
eled is only about 10 Mpc, or 5% of the size of the sim- 
ulation box, and the largest distance traveled is about 15% 
of this size. So we may assume that in the optimal assign- 
ment distances between paired positions will be limited. We 
define then a critical distance dcrit and specify that a final 
position Xi and an initial position form a valid pair only 
if they are within less than dcrit from each other. This crit- 
ical distance must be adjusted carefully: if it is too small, 
we risk excluding the optimal assignment; if it is taken too 
large, the benefit of the sparse description is lost. 

However, the saving in computing time achieved by 
sparse description has to be paid for in storage space: to 
store the set A of valid pairs, storage of size |.4| — fN"^ is 
needed, which takes us back to the 0{N^) storage require- 
ment. We have explored two solutions to this problem. 

1. Use a dense description nevertheless, i.e. the one 
where all pairs {i,j) are valid and there is no need to store 
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the set A. The auction program is easily adapted to this 
case (in fact this simplifies the code). However, we forfeit 
the saving in time provided by the sparse structure. 

2. The sparse description can be preserved if the set 
of valid pairs is computed on demand rather than stored. 
This is easy if initial positions fill a uniform cubic grid, the 
simplest discrete approximation to the initial quasi-uniform 
distribution of matter in the reconstruction problem. Thus, 
for a given final position Xi, the valid pairs correspond to 
points of the cubic lattice that lie inside a sphere of radius 
dcrit centered at Xi, so their list can be generated at run 
time. 

Fig. El gives the computing time as a function of the 
number of points N used in the assignment problem. Shown 
are the dense and sparse versions of the auction algorithm 
(in the latter, t he critical distance square d was taken equal 
to 200) and the lBurkard fc Derigj lll980l) algorithm, which 
ranked the next fastest in our experiments. The N initial and 
final positions are chosen from the file generated by an A''- 
body simulation described in Section|21 the choice is random 
except for the sparse algorithm, in which the initial positions 
are required to fill a cubic lattice. Hence, the performance 
of the sparse auction algorithm shown in the figure is not 
completely comparable to that of the two other algorithms. 

It is evident that the differenc e in computing tim e be- 
tween the dense auction and the iBurkard fc Derie^ algo- 
rithms steadily increases. In the vicinity of A'^ = 10^, the 
dense auction algorithm is about 10 times faster than the 
other one. For the sparse version, the decrease in computing 
time is spectacular: as could be expected, the ratio of com- 
puting times for the two versions of the auction algorithm 
is of the order of /. For large A'^, the 0{N'^) asymptotic of 
the computing time is quite clear for the sparse auction al- 
gorithm. For two other algorithms, similar asymptotic was 
found for larger A'^ in other experiments (not shown). 

In all three cases shown, the initial positions fill a con- 
stant volume while A' is varied. This is what we call constant- 
volume computations. In the sparse case, this results in a 
constant filling factor, equal to the ratio of the volume of 
the sphere with radius dcrit to the volume occupied by the 
initial positions. Here this filling factor is about / = 0.019. 
Another choice, not shown in the figure, is that of constant- 
density computations, when the initial positions are taken 
from a volume whose size increases with A'^. In this case the 
time dependence of algorithms for large A'^ is of the order 
of N^-^. 

We finally observe that the sparse auction algorithm ap- 
plied to the MAK reconstruction requires 5 hours of single- 
processor CPU time on a 667 MHz COMPAQ/DEC Alpha 
machine for 216,000 points. 



5 TESTING THE MAK RECONSTRUCTION 

In this section we present results of our testing the MAK 
reconstruction against data of cosmological A^-body simu- 
lations. In a typical simulation of this kind, the dark mat- 
ter distribution is approximated by A' particles of identical 
mass. Initially the particles are put on a uniform cubic grid 
and given velocities that form a realization of the primor- 
dial velocity field whose statistics is prescribed by a cer- 
tain cosmological model. Trajectories of particles are then 



Figure 7. Af-body simulation output in the Eulerian space used 
for testing our reconstruction method (shown is a projection 
onto the x-y plane of a 10% slice of the simulation box of size 
200h-l Mpc). Points are highlighted in yellow when reconstruc- 
tion fails by more than 6.25 Mpc, which happens mostly in 
high-density regions. 

computed according to the Newtonian dynamics in a co- 
moving frame, using periodic boundary conditions. The re- 
construction problem is therefore to recover the pairing be- 
tween the initial (Lagrangian) positions of the particles and 
their present (Eulerian) positions in the A'-body simulation, 
knowing only the set of computed Eulerian positions in the 
physical space. 

We test our reconstruction against a simulation of 
128"^ particles in a box of 200 Mpc size (where h 
is the Hubble parameter in units of 100 km s~^ Mpc~^) 
performed using the adaptive P ^M code HYDRA 
JCouchman. Thomas fc Pearcd ITooil) ." A ACDM cosmo- 
logical model is used with parameters Q,m = 0.3, Q.a = 0.7, 
h = 0.65, (T8 = 0.9.^^ The value of these parameters within 
the model are determined by fitting the observed cosmic mi- 
crowave background (CMB) spectrum. The output of the 
TV-body simulation is illustrated in Fig. by a projection 
onto the x-y plane of a 10% slice of the simulation box. 

Since the simulation assumes periodic boundary condi- 
tions, some Eulerian positions situated near boundaries may 
have their Lagrangian antecedents at the opposite side of the 

In a flavour of A'-body codes called particle-mesh (PM) codes, 
Newtonian forces acting on particles are interpolated from the 
gravitational field computed on a uniform mesh. In very dense 
regions, precision is increased by adaptively refining the mesh and 
by direct calculation of local particle-particle (PP) interactions; 
codes of this type are correspondingly called adaptive P^M. 

The use of a ACDM model instead of the model without a 
cosmological constant (Appendix leads to some modifications 
in basic equations but does not change formulas used for the MAK 
reconstruction. 

Data of the first year Wil kinson M icrowave Anisotropy Probe 
iSoergel et al.lbOO.'j: see also iBridle et al. 2003) suggest a value 
(Tg = 0.84 it 0.04, marginally smaller than the one used here. This 
may slightly extend the range of scales favourable for the MAK 
reconstruction. 
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simulation box. Suppressing the resulting spurious large dis- 
placements is crucial for successful reconstruction. Indeed, 
for a typical particle displacement of 1/20 the box size, spu- 
rious box- wide leaps of 1% of the particles will generate a 
contribution to the quadratic cost I34II four times larger than 
that of the rest. To suppress such leaps, for each Eulerian 
position that has its antecedent Lagrangian position at the 
other side of the simulation box, we add or subtract the box 
size from coordinates of the latter (in other words, we are 
considering the distance on a torus) . In what follows we refer 
to this procedure as the periodicity correction. 

We first present reconstructions for three samples of 
particles initially situated on Lagrangian subgrids with 
meshes given by Ax = 6.25/i"^Mpc, Ax/2 and Ax/A. To 
further reduce possible effects of the unphysical periodic 
boundary condition, we truncate the data by discarding 
those points whose Eulerian positions are not within the 
sphere of radius 16Ax placed at the centre of the simulation 
box (for the largest Ax its diameter coincides with the box 
size). The problem is then confined to finding the pairing be- 
tween the remaining Eulerian positions and the set of their 
periodicity-corrected Lagrangian antecedents in the A'^-body 
simulation. 

The results are shown in Figs. ISHTTl The main plots 
show the scatter of reconstructed vs. simulation Lagrangian 
positions for the same Eulerian positions. For these diagrams 
we introduce a 'quasi-periodic projection' 

q=iqi + V2q2 + ^93)7(1 + V2 + VS) (42) 

of the vector q, which ensures a one-to-one correspondence 
between g- values and points on the regular Lagrangian grid. 
The insets are histograms (by percentage) of distances, in 
reconstruction mesh units, between the reconstructed and 
simulation Lagrangian positions; the first darker bin, slightly 
less than one mesh in width, corresponds to perfect recon- 
struction (thereby allowing a good determination of the pe- 
culiar velocities of galaxies). 

With the mesh size Ax, Lagrangian positions of 62% of 
the sample of 17,178 points are reconstructed perfectly and 
about 75% are placed within not more than one mesh. With 
the Ax/2 grid, we still have 35% of exact reconstruction out 
of 19,187 points, but only 14% for the Ax/A grid with 23,111 
points. 

We also performed a reconstruction on a random sample 
of 100,000 Eulerian positions taken with their periodicity- 
corrected Lagrangian antecedents out of the whole set of 
128'^ particles, without any restrictions. This reconstruction, 
with the effective mesh size (average distance between neigh- 
bouring points) of 4.35^~^Mpc, gives 51% of perfect recon- 
struction fFig. IIH . 

We compared these results with those of the 
PIZA reconstruction method (see Section 14.11 and 
ICroft &: Gaztaiiagall99'^ , which gives a 2-monotone but not 
necessarily optimal pairing between Lagrangian and Eule- 
rian positions. We applied the PIZA method on the Aa; grid 
and obtained typically 30-40% exactly reconstructed posi- 
tions, but severe non-uniqueness: for two different seeds of 
the random generator used to set up the initial tentative as- 
signment, only about half of the exactly re constructed posi- 
tions were the same (see figs. 3 and 7 of iMohavaee et all 
l|2223) for an illustration). We also implemented a mod- 
ification of the PIZA method establishing 3-monotonicity 
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simulation coordinate 

Figure 8. Test of the MAK reconstruction for a sample of 
N' = 17, 178 points initially situated on a cubic grid with mesh 
Ax = 6.25 Mpc. The scatter diagram plots true versus recon- 
structed initial positions using a quasi-periodic projection which 
ensures one-to-one correspondence with points on the cubic grid. 
The histogram inset gives the distribution (in percentages) of dis- 
tances between true and reconstructed initial positions; the hori- 
zontal unit is the sample mesh. The width of the first bin is less 
than unity to ensure that only exactly reconstructed points fall 
in it. Note that more than sixty percent of the points are exactly 
reconstructed. 
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Figure 9. Same as Fig. 13 but with N' = 19,187 and a sample 
mosh of Ax/2 = 3.125 Mpc. Exact reconstruction is down to 
35%. 

(monotonicity with respect to interchanges of 3 points in- 
stead of pairs) and checked that it does not give a significant 
improvement over the original PIZA. 

In comoving coordinates, the typical displacement of 
a mass element is about 1/20 the box size, that is about 
Wh~^ Mpc. This is not much larger than the coarsest grid 
of 6.25 /i"^ Mpc used in testing MAK which gave 62% of 
exact reconstruction. Nevertheless there are 18 other grid 
points within lOft^^Mpc of any given grid point, so that 
this high percentage cannot be trivially explained by the 
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Figure 10. Same as Fig.l^Jbut with N' = 23, 111 and a sample 
mesh of Ax/4 = 1.55 Mpc. Exact reconstruction is down to 
14%. 
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Figure 11. Same as Fig. Isl with TV' = 10^ points selected at 
random, neighbouring points being typically 4.35 Mpc apart. 
Exact reconstruction is in excess of 50%. 



smallness of the displacement. Note that without the peri- 
odicity correction, the percentage of exact reconstruction for 
the coarsest grid degraded significantly (from 62% to 45%) 
and the resulting cost was far from the true minimum. 

For real catalogues, reconstruction has to be performed 
for galaxies whose positions are specified in the redshift 
space, where they appear to be displaced radially (along the 
line of sight) by an amount proportional to the radial com- 
ponent of the peculiar velocity. Thus, at the present epoch, 
the redshift position s of a mass element situated at the 
point X in the physical space is given by 



s = X + xP (v ■ x) 



(43) 



where v is the peculiar velocity in the comoving coordi- 
nates X and the linear growth factor time r, x denotes the 
unit normal in the direction of x , and the parameter j3 equals 
0.486 in our ACD M model. 

Following IValentine et alJ ll2000l : see also 




Figure 12. Test of the redshift-space variant of the MAK recon- 
struction based on the same data as Fig. |S] The circular redshift 
map (violet points) corresponds to the same physical-space slice 
as displayed in Fig. |7| (the observer is taken at the center of the 
simulation box). Points are highlighted in red when reconstruc- 
tion fails by more than one mesh. 



iMonaco fc EfstathiotJ Il99flh . we use the Zel'dovich ap- 
proximation (ZA) to render our MAK quadratric cost 
function in the s variable. As follows from pi^ . in this 
approximation the peculiar velocity is given by 



V = —{x — q). 

T 



At the present time, since tq 
gives 



(44) 

1, this together with l|43^ 



{s — q) ■ X = {1 + f3){x — q) ■ X, 



+ P{P + 2){{x-q)-xy 



(45) 
(46) 



Combining now these two equations and using the fact that, 
by 14311 . the vectors x and s are coUinear and therefore x = 
±s, we may write the quadratic cost function as 



2(/3 + l)2^^ 



(47) 



The redshift-space reconstruction is then in principle re- 
duced to the physical-space reconstruction. Note however 
that the redshift transformation of Eulerian positions may 
fail to be one-to-one if the peculiar component of velocity 
field in the proper space coordinates exceeds the Hubble ex- 
pansion component. This undermines the simple reduction 
outlined above for catalogues confined to small distances. 

We have performed a MAK reconstruction with the 
redshift-modified cost function 1471 . The redshift positions 
were computed for the simulation data with peculiar veloc- 
ities smoothed over a sphere with radius of 1/100 the box 
size {2h~^ Mpc). This reconstruction led to 43% of exactly 
reconstructed positions and 60% which are within not more 
than one Ax mesh from their correct positions (see Fig. 1121 
a scatter diagram is omitted because it is quite similar to 
that in Fig. I^J. A comparison of the redshift-space MAK 
reconstruction with the physical-space MAK reconstruction 
shows that almost 50% of exactly reconstructed positions 
correspond to the same points. This test shows that the 
MAK method is robust with respect to systematic errors 
introduced by the redshift transformation. 

Our results demonstrate the essentially potential char- 
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acter of the Lagrangian map above ~ 6 Mpc (within the 
ACDM model) and perhaps at somewhat smaller scales. 

Although it is not our intention in this paper to actually 
implement the MAK reconstruction on real catalogues, a few 
remarks are in order. The effect of the catalogue selection 
function can be handled by standard techniques; for instance 
one can assign each galaxy a 'm ass' inversely proportional 
to the catalog selection function iNusser fc Branchinilliooot 
IValentine et alJl200nl : iBranchini et al.ll2W .^- Biasing can be 
taken into account in a similar manner ijNusscr & Branchini 
I2OO0I) . Both these modifications and the natural scatter of 
masses in the observational catalogues require that mas- 
sive objects be represented by clusters of multiple Eule- 
rian points of unit mass (with the correspondingly increased 
number of points on a finer grid in the Lagrangian space), 
which reduces the problem to a variant of the usual assign- 
ment. We also observe that real catalogues involve trun- 
cation, that is data available only over a finite region. As 
already discussed in Section [3.41 this is not a serious prob- 
lem provided a sufficiently large patch is available. Actu- 
ally, as noted earlier in this Section, the data used in testing 
have been truncated spherically, without significantly affect- 
ing the quality of the reconstruction. 

In the redshift-space modification, more accurate de- 
termination of peculiar velocities can be done using second- 
order Lagrangian perturbation theory. Note also that, for 
the observational catalogues, the motio n of the local group 
itself should also be accounted for ijTavlor fc Valentina 

Iliii). 



6 RECONSTRUCTION OF THE FULL 
SELF-GRAVITATING DYNAMICS 

The MAK reconstruction discussed in Sections |3 and |1| was 
performed under the assumption of a potential Lagrangian 
map and of the absence of multi-streaming. The tests done 
in Section 1^ indicate that potentiality works well at scales 
above 6/i~^Mpc, whereas multi-streaming is mostly be- 
lieved to be unimportant above a few megaparsecs. There 
could thus remain a substantial range of scales over which 
the quality of the reconstruction can be improved by re- 
laxing the potentiality assumption and using the full self- 
gravitating dynamics. Here we show that, as long as the dy- 
namics can be described by a solution to the Euler-Poisson 
equations, the prescription of the present density field still 
determines a unique solution to the full reconstruction prob- 
lem. We give only the main ideas, technical details being 
left for Ap pendix [p] (a, m athematically rigorous proof may 
be found in lLoeoerl l2003h ). In order to make the exposition 
self-contained, we also give in Appendix^lan elementary in- 
troduction to convexity and duality which are used for the 
derivation (and also elsewhere in this paper). 

We shall start from an Eulerian variational formulation 
of the Euler-Poisson equations in an Einstein-de Sitter uni- 
ve rse, which is an adaptat ion of a variational principle given 
bv lGiavalisco et al.l l)l993h . We minimize the action 

/=^^°dry d^a3r^/^(p|«|2 + ^|V.v'g|'), (48) 

under the following four constraints: the Poisson equation 
©, the mass conservation equation (|5J and the boundary 



conditions that the density field be unity at r = and pre- 
scribed at the present time t — tq. The constraints can 
be handled by the standard method of Lagrange multipli- 
ers (here functions of space and time) , which allows to vary 
independently the fields p, (fg and v. The vanishing of the 
variation in v gives v = t~^^'^V ^d, where 9{x,t) is the 
Lagrange multiplier for the mass conservation constraint. 
Hence, the velocity is curl-free. The vanishing of the varia- 
tion in p gives then 

^^^ + ^|V=^^l' + |:^-0- (49) 

By taking the gradient, this equation goes over into the mo- 
mentum equation Q, repeated here for convenience: 

3 

drV + {v-V^)v = - — {v + V^^s). (50) 

AT 

It is noteworthy that, if in the action we replace 3/2 
both in the exponent of r and in the gravitational energy 
term by 3q/2, we obtain l|50^ but also with a 3a/(2r) factor 
in the right-hand side. The Zel'dovich approximation and 
the associated MAK reconstruction amount clearly to set- 
ting Q = 0, so as to recover the 'free-streaming action' 

I dr J d'x p\v\\ (51) 

whose minimization is easily shown to be equivalent to that 
of the quadratic cost function II23II . 

Assuming the action l|48|l to be finite, existence of a min- 
imum is mostly a consequence of the action being manifest- 
edly non-negative. Here it is interesting to observe that the 
Lagrangian, which is the difference between the kinetic en- 
ergy and the potential energy, is positive whereas the Hamil- 
tonian which is their sum does not have a definite sign. As 
a consequence, our two-point boundary problem is, as we 
shall see, well posed but the initial- value problem for the 
Euler-Poisson system is not well posed since formation of 
caustics after a finite time cannot be ruled out.^* 

Does the variational formulation imply uniqueness of 
the solution? This would be the case if the action were a 
strictly convex functional (see Appendix IClfl . which is guar- 
anteed to have one and only one minimum. The action as 
written in 1481 is not convex in the p and v variables, but 
can be rendered so by introducing the mass fiux J = pv; 
the kinetic energy term becomes then |Jp/(2p), which is 
convex in the J and p variables. 

Strict convexity is particularly cumbersome to estab- 
lish, but there is an alternative way, known as duality: by 
a Legendre-like transformation the variational problem is 
carried into a dual problem written in terms of dual vari- 
ables; the minimum value for the original problem is the 
maximum for the dual problem. It turns out that the dif- 
ference of these equal values can be rewritten as a sum of 
non-negative terms, each of which must thus vanish. This is 
then used to prove (i) that the difference between any two 
solutions to the variational problem vanishes and (ii) that 
any curl-free solution to the Euler-Poisson equations with 
the prescribed boundary conditions for the density also min- 



If wc had considered electrostatic repulsive interactions the 
conclusions would be reversed. 
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imizes the action. All this together establishes uniqueness. 
For details see Appendix IHl 

Several of the issues raised in connection with the 
MAK reconstruction appear in almost the same form for the 
Euler-Poisson reconstruction. First, we are faced again with 
the problem that, when reconstructing from a finite patch of 
the present universe, we need either to know the shape of the 
initial domain or to make some hypothesis as to the present 
distribution of matter outside this patch. Second, just as for 
the MAK reconstruction, the proof of uniqueness still holds 
when the present density poix) has a singular part, that is, 
when some matter is concentrated. Again, we shall have full 
information on the initial shape of collapsed regions but not 
on the initial fiuctuations inside them. The particular solu- 
tion obtained from the variational formulation is the only 
solution which stays smooth for all times prior to tq. 

We also note that, at this moment and probably for 
quite some time, 3D catalogues sufficiently dense to allow 
reconstruction will be limited to fairly small redshifts. Even- 
tually, it will however become of interest to perform recon- 
struction 'along our past light-cone' with data not all at 
Tq. The variational approach can in principle be adapted to 
handle such reconstruction. 

In previous sections we have seen how to implement re- 
construction using MAK, which is equivalent to using the 
simplified action 15111 . Implementation using the full Euler- 
Poisson action 1481 is mostly beyond the scope of this paper, 
but we shall indicate some possible directions. In principle 
it should be possible to adapt to the Euler-Poisson recon- 
struction the method of the augmented Lagrangian which 
has been applied to the two-dimen sional Monge-Ampere 
equation jBenamou fc Bren ier"200Cf) . An alternative strat- 
egy, which allows reduction to MAK- type problems, uses the 
idea of 'kicked burgulence' l(Bec. Frisch fc Khaninl l2000i in 
which, in order to solve the one or multi-dimensional Burg- 
ers equation 

drV + {v ■Vx)v = f(x,T), V = -Vx<Pw, (52) 

one approximates the force by a sum of delta-functions in 
time: 

f{x,T)^Y.S{T-n)g^{x). (53) 

i 

In the present case, the g^ix) are proportional to the right- 
hand side of 15011 evaluated at the kicking times r^. The ac- 
tion becomes then a sum of free-streaming Zel'dovich-type 
actions plus discrete gravitational contributions stemming 
from the kicking times. Between kicks one can use our MAK 
solution. At kicking times the velocity undergoes a discon- 
tinuous change which is related to the gravitational poten- 
tial (and thus to the density) at those times. The densities 
at kicking times can be determined by an iterative proce- 
dure. The kicking strategy also allows to do redshift-space 
reconstruction by applying the redshift-space modified cost 
(Section at the last kick. 



7 COMPARISON WITH OTHER 
RECONSTRUCTION METHODS 

Reconstruction started with Peebles' (1989) work, in which 
he compared reconstructed and measured peculiar veloci- 



ties for a small number of Local Group galaxies, situated 
within a few Mpc. The focus of reconstruction work has 
now moved t o tackling the rapidly gr owing large 3D sur- 
veys (see, e.g. iFrieman fc Szalavll200Cl) . It is not our inten- 
tion here to review all the work on reconstruction;^^ rather 
we shall discuss how some of the previously used meth- 
ods can be reinterpreted in the light of the optimization 
approach to reconstruction. For convenience we shall di- 
vide methods into perturbative (Section 17.1^ . probabilistic 
(Section 17.21 1. and vsjiational (S ection 17.311 . Methods such 
as POTENT jDekel et al.lll99Ci) . whose purpose is to ob- 
tain the full peculiar velocity field from its radial compo- 
nents using the (Eulerian) curl-free property, are not directly 
within our scope. Note that in its original Lagrangian form 
( Bcrtschinecr fc Dckcl 1989; Dokcl ct al. 1990) POTENT 
was assuming a curl-free velocity in Lagrangian coordinates, 
an assumption closely related to the potential assumption 
made for MAK, as already pointed out in Section f3.1l Even 
closer is the relation between MAK and the PIZA method of 
Croft fc Gaztanaaa (1997), discussed in Section [7.31 which 
is also based on minimization of quadratic action. 

7.1 Perturbative methods 

iNusser fc DekeJ lll992h have proposed using the Zel'dovich 
approximation backwards in time to obtain the initial ve- 
locity fiuctuations and thus (by slaving) the density fiuc- 
tuations. Schematically, their procedure involves two steps: 
(i) obtaining the present potential velocity field and (ii) in- 
tegrating the Zel'dovich-Bernouilli equation back in time. 
Using the equality (in our notation) of the velocity and grav- 
itational potentials, they point out that the velocity poten- 
tial can be computed from the present density fluctuation 
field by solving the Poisson equation. This is a perturba- 
tive approximation to reconstruction in so far as it replaces 
the Monge-Ampere equation H19|l by a linearized form. In- 
deed, when using the Zel'dovich approximation we have 
q — X — TV = x + TVxfvix). We know that q = Va;©(a;) 
with B satisfying the Monge-Ampere equation. The latter 
can thus be rewritten as 

det [S^J + rV:,,V^j((3v(a;)) = p{x), (54) 

where 5ij denotes the identity matrix. If we now use the 
relation det{5ij + eAij) = 1 -|- e ^ . Au -\- 0{e^) and truncate 
the expansion at order e, we obtain the Poisson equation 

TVl,ip^{x) ^ p{x) - 1 = 5{x). (55) 

Of course, in one dimension no approximation is needed. 
From a physical point of view, equating the velocity and 
gravitational potentials at the present epoch amounts to us- 
ing the Zel'dovich approximation in reverse and is actually 
inconsistent with the forward Zel'dovich approximation: the 
slaving which makes the two potentials equal initially does 
not hold in this approximation at later epochs. Replacing 
the Monge-Ampere equation by the Poisson equation is not 
consistent with a uniform initial distribution of matter and 
will in general lead to spurious multi-streaming in the ini- 
tial distribution. Of course, if the present-epoch velocity field 

For a comparison of six different techniques, see 
iNaravanan fc Crofj <199Sl) . 
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happens to be known one can try applying the Zel'dovich 
approximation in reverse. Nusser and Dekel observe that cal- 
culating the inverse Lagrangian map hy q = x — tv does not 
work well (spurious multi-streaming appears) and instead 
integrate back in time the Zel'dovich-Bernouilli equation^" 

dt<fiv = ^ C^xipvf , (56) 

which is obviously equivalent to the Burgers equation Ijli-ifl 
with the viscosity 1/ — 0. One way of performing this re- 
verse integration, which guarantees the absence of multi- 
streaming, is to use the Legendre transformation l|18|l to 
calculate ^{q) from Q{x) = \x\'^ /2 — T(pv{x) and then ob- 
tain the reconstructed initial velocity field as 

«in(<7) = «o(V,<l>(q)). (57) 

This procedure can however lead to spurious shocks in the 
reconstructed initial conditions, due to inaccuracies in the 
present-epoch velocity data, unless the data are suitably 
smoothed . Final ly, the improved reconstruction method of 
GramannI lll993l) can be viewed as an approximation to 
the Monge-Ampere equation beyond the Poisson equation 
which captures part of the nonlinearity. 

7.2 Probabilistic methods 

IWeinberel l)l992h presents an original approach to recon- 
struction, which turns out to have hidden connections to 
optimal mass transportation. The key observations in his 
'Gaussianization' technique are the following: (i) the initial 
density fiuctuations are assumed to be Gaussian, (ii) the 
rank order of density values is hardly changed between initial 
and present states, (iii) the bulk displacement of large-scale 
features during dynamical evolution can be neglected. As- 
sumption (i) is part of the standard cosmological paradigm. 
Assumption (iii) can of course be tested in A^-body sim- 
ulations. As we have seen in Section |3 a displacement of 
10^~^Mpc is typical and can indeed be considered small 
compared to the size of the simulation boxes (64/i-iMpc m 
Weinberg's simulations and 200/i~^Mpc in ours). Assump- 
tion (ii) means that the correspondence between initial and 
present values of the density p (or of the contrast S = p—1) 
is monotone. This map, which can be determined from the 
empirical present data, can then be applied to all the data to 
produce a reconstructed initial density field. Finally, by run- 
ning an A^-body simulation initialized on the reconstructed 
field one can test the validity of the procedure, which turns 
out to be quit e good and can be improved further by hy- 
brid methods ijNaravanan fc Weinberg! \l99fi iKolatt et all 
Il996^ comb ining G aussianization wit h the perturbativ e ap- 
proaches of iNusser fc Dokcl (1991) or lGramannI ijim^t) . 

This technique is actually connected wi th mas s trans- 
portatio n: starting wit h the work of Frechet Jl957atll957bt 
see also lRachevlll983^ . probabilists have been asking the 
following question: given two random variables mi and m2 
with two laws, say PDFs pi and p2, can one find a joint 
distribution of (mi,m2) with PDF pi2(mi,m2) having the 

^'^ In the non-cosmological literature this equation is usually 
called Hamilton-Jacobi i n the context of analytical mecha nics 
jLandau &: Lifshit j and Kardar—Parisi— Zhang ll98fift in 

condensed matter physics. 



following properties: (i) pi and p2 are the marginals, i.e. 
when P12 is integrated over m2 (respectively, mi) one re- 
covers pi (respectively, P2), (ii) the correlation {mim2) is 
maximum? Since (mf) and (m2) are obviously prescribed 
by the constraint that we know pi and P2, maximizing the 
correlation is the same as minimizing the quadratic distance 
{(mi — m2)'^) . This is precisely an instance of the mass trans- 
portation problem with quadratic cost, as we defined it in 
Section 13.31 As we know, the optimal solution is obtained 
by a map from the space of mi values to that of m2 val- 
ues which is the gradient of a convex function. If mi and 
m2 are scalar variables, the map is just monotone, as in the 
Gaussianization method (in the discrete setting this was al- 
ready observed in Section 14. H . Hence Weinberg's method 
may be viewed as requiring maximum correlation (or mini- 
mum quadratic distance in the above sense) between initial 
and present distributions of density fiuctuations. 

In principle the Gaussianization method can be ex- 
tended to multipoint distributions, leading to a difficult mul- 
tidimensional mass transportation problem which can be 
discretized into an assignment problem just as in Section im 
The contact of the maximum correlation assumption to the 
true dynamics is probably too flimsy to justify using such 
heavy machinery. 

T.3 Variational methods 

All vari ational a pproa ches to reconstruction, starting with 
that of |Peeble3 lll989f) . have common features: one uses a 
suitable Lagrangian and poses a two-point variational prob- 
lem with boundary conditions prescribed at the present 
epoch by the observed density field, and at early times by 
requiring a quasi-uniform distribution of matter (more pre- 
cisely, as we have seen in Section f2. II by requiring that the 
solutions not be singular as r — > 0). 

The Path Interchange Zel'dovich Approximation 
(PIZA) method of lGroft fc Gaztanaeal Jl997l) and our MAK 
reconstruction techniques use a free-streaming Lagrangian 
in linear growth rate time. As we have seen in Sec- 
tion 13.11 this amounts to assuming adhesion dynamics. 
Once discretized for numerical purposes, the variational 
problem becomes an inst ance of the assignment problem. 
ICroft fc Gaztanagal il997ll have proposed a restricted pro- 
cedure for solving it, which does not account for the La- 
grangian potentiality and yields non-unique approximate 
solutions. As we have seen in Sections 2] and |K1 the exact 
and unique solution can be found with reasonable CPU re- 
sources. 

Turning now to the Peebles least action method, let 
us first describe it schematically, using our notation. In its 
original formulation it is applied to a discrete set of galaxies 
(assumed of course to trace mass) in an Einstein-de Sitter 
universe. The action, in our notation, can be written as 

r ^ ^ ( "liT^ \dXi 1^ 

I = dr^l^^^l — I 

where m^ is the mass and the comoving coordinate of ith 
galaxy (see also iNusser fc Branchinill200iOf) . This is supple- 
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Figure 13. A schematic demonstration of Peebles' reconstruc- 
tion of the trajectories of the members of the local neighbourhood 
using a variational approach based on the minimization of Euler— 
Lagrange action. The arrows go back in time, starting from the 
present and pointing towards the initial positions of the sources. 
In most cases there is more than one allowed trajectory due to 
orbit crossing (closely related to the multi-streaming of the under- 
lying dark matter fluid). The pink (darker) orbits correspond to 
taking the minimum of the action whereas the yellow (brighter) 
orbits were obtained by taking the saddle-point solution. Of par- 
ticular interest is the orbit of N6822 which in the former solution 
is on its first approach towards us and in the second solution is in 
its passing orbit. A better agreement between the evaluated and 
observed velocities was shown to correspond to the saddle-point 
solution. 



mented by the boundary condition that the present positions 
of the galaxies are known and that the early-time velocities 
satisfy^^ 



3/2 



dxi 







for 



0. 



(59) 



This particle approach was extended by iGiavalisco et alJ 
l)l993h to a continuous distribution in Eulerian coordinates 
and leads then to the action analogous to 14811 which we have 
used in Section |S| The procedure also involves a 'Galerkin 
truncation' of the particle trajectories to finite sums of trial 
functions of the form 



<W = <{r,) + Y.ClJ^{T), 

n=0 

/n(r) = t"{to-t), n = 0,l,...,7V- 1. 



(60) 
(61) 



The reconstructed peculiar velocities for the Local Group 
were used by Peebles to calibrate the Hubble and den- 
sity parameters, which turned out to differ from the previ- 
ously assumed values. However the peculiar velocity of one 
dwarf galaxy, N68 22, failed to ma tch the observed value (see 
Fig. ESJ. This led IPeeblesI Jl990l) to partially relax the as- 
sumption of minimum action, allowing also for saddle points 



This condition, which is written o?dxi/dt — > in Peebles' 
notation, ensures the vanishing of the corresponding boundary 
term after an integration by parts in the time variable. 



in the action. Somewhat better agreement with observations 
is then obtained, but at the expense of lack of uniqueness. 

In the context of the present approach, various remarks 
can be made. The boundary condition l|59^ is trivially sat- 
isfied if the velocities dx/dr remain bounded. Actually, we 
have seen in Section l2. II that, as a consequence of slaving, 
the velocity has a regular expansion in powers of r, which 
implies its boundedness as r ^ 0. The important point is 
that the function /n(T) appearing in H6()|l should be expand- 
able in powers of r, as is the case with the ansatz H61^ . 

In Section]^ we have established uniqueness of the re- 
construction with a prescribed present density and under the 
assumption of absence of multi-streaming (but we allow for 
mass concentrations). This restriction is meaningful only in 
the continuous case: in the discrete case, unless the particles 
are rather closely packed, the concept of multi-streaming is 
not clear but there have been attemp ts to relate uniqueness 
to absence of ' orbit crossing' (see, e.g. . IGiavalisco et alll993 : 
IWhitin J200ffl . Of course, at the level of the underlying dark 
matter, multi-streaming is certainly not ruled out at suffi- 
ciently small scales; at such scales unique reconstruction is 
not possible. 

In the truly discrete case, e.g. when considering a dwarf 
galaxy, there is no reason to prefer the true minimum action 
solution over any other stationary action solution. 



8 CONCLUSION 

The main theoretical result of this paper is that reconstruc- 
tion of the past dynamical history of the Universe, know- 
ing only the present spatial distribution of mass, is a well- 
posed problem with a unique solution. More precisely, re- 
construction is uniquely defined down to those scales, a 
few megaparsecs, where multi-streaming becomes impor- 
tant. The presence of concentrated mass in the form of 
clusters, filaments, etc is not an obstacle to a unique dis- 
placement reconstruction; the mass within each such struc- 
ture originates from a collapsed region of known shape but 
with unknown initial density and velocity fiuctuations in- 
side. There are of course practical limitations to reconstruc- 
tion stemming from the knowledge of the present mass dis- 
tribution over only a limited patch of the Universe; these 
were discussed in Section f3. 41 

In this paper we have also presented in detail and tested 
a reconstruction method called MAK which reduces recon- 
struction to an assignment problem with quadratic cost, for 
which effective algorithms are available. MAK, which is ex- 
act for dynamics governed by the adhesion model, works 
very well above 6 Mpc and can in principle be adapted 
to full Euler-Poisson reconstruction. 

We note that a very common method for testing ideas 
about the early Universe is to take some model of early den- 
sity fluctuations and then run an A''-body simulations with 
assumed cosmological parameters until the present epoch. 
Confrontation with the observed statistical properties of the 
present Universe helps then in selecting plausible models and 
in narrowing the choice of cosmological parameters. This for- 
ward method is conceptually very different from reconstruc- 
tion; the latter not only works backward but, more impor- 
tantly, it is a deterministic method which gives us a detailed 
map of the early Universe and how it relates to the present 
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one. Reconstruction thus allows us to obtain the peculiar 
velocities of galaxies and is probably the only method which 
can hope to do this for a large number of galaxies. In those 
instances were we have partial information on peculiar veloc- 
ities (from independent distance meas urements), e .g. for the 
Near By Galaxies (NBG) catalogue of lTuUvl l)l988|) . such in- 
formation can be used to calibrate cosmological parameters 
or to provide additional constraints, which are in principle 
redundant but can improve the quality. 

The detailed reconstruction of early density fluctua- 
tions, which will become possibl e using large 3D survey s 
such as 2dF and SDSS (see, e.g.. iFrieman fc Szalavlf200(t . 
will allow us to test such assumptions as the Gaussianity of 
density fluctuations at decoupling. Note however that such 
reconstruction gives us full access only to the complement of 
collapsed regions; any statistical information thus obtained 
will be biased, roughly by overemphasizing underdense re- 
gions. 

Finally we have no reason to hide the pleasure we ex- 
perience in seeing this heavenly problem bring together and 
indeed depend crucially on so many different areas of math- 
ematics and physics, from fluid dynamics to Monge- Ampere 
equations, mass transportation, convex geometry and com- 
binatorial optimization. Probably this is the first time that 
one tackles the three-dimensional Monge-Ampere equation 
numerically for practical purposes. As usual, we can ex- 
pect that the techniques, here applied to cosmic reconstruc- 
tion, will find many applications, for example to the optimal 
matching of two holographic or tomographic images or to the 
correction of images in multi-dimensional colour space. 
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APPENDIX A: EQUATIONS OF MOTION IN 
AN EXPANDING UNIVERSE 

On distances covered by present and forthcoming red- 
shift galaxy catalogues, the Newtonian description consti- 
tutes a realistic approximation to the dynamics of self- 
gravi t ating cold dark matter filling the Universe (jPeeblej 
Il980l: IColes fc LucchinI l2002h . This description gives, in 
proper space coordinates denoted here by r and cosmic 
time t, the familiar Euler-Poisson system for the den- 
sity Q{r,t), velocity U{r,t) and the gravitational poten- 
tial (/)(r, t): 



dtU+{U-\7r)U 

dtg + Vr- (gU) 



= 0, 

= inGg, 



(A.l) 
(A.2) 
(A.3) 



where G is the gravitation constant. 

In a homogeneous isotropic universe, the density and 
velocity fields take the form 



gir,t)^g{t), U{r,t) 



H(t)r = -^r. 



(A.4) 



Here the coefficient H{t) is the Hubble parameter, and a{t) is 
the expansion scale factor defined so that integration of the 
velocity field r — U{r,t) = H{t)r yields r — a{t)x, where 
X is called the comoving coordinate. 

The background density g{t) gives rise to the back- 
ground gravitational potential 0g, which by ijA.lfl and l|A.4^ 
satisfies 



a 

— —r. 

a 



(A.5) 



For the background density, mass conservation IIA.2t gives 
then 



ga^ = £10, 



(A.6) 

where ~ gito) with to the present epoch and a{to) is 
normalized to unity. Eqs. IIA.5L IIA.6L and <A.3ll imply the 
Friedmann equation for a{t): 

4 1 

a = --7tG£>o^ (A.7) 
with conditions posed at t = to: 

a{to) = 1, d{to) ^Ho>0, (A.8) 

where Ho is the present value of the Hubble parameter, pos- 
itive for an expanding universe. 

For simplicity we restrict ourselves to the case of the 
critical density, corresponding to the flat, matter-dominated 
Einstein-de Sitter universe (without a cosmological con- 
stant): 



Qo = 



3Hl 

8nG 



(A.9) 



and adjust the origin of the time axis such that the solution 
takes the form of a power law 

a{t) = (Ij (A.IO) 

with Ho = 2/(3to) and = l/{6nGtl). 

The observed Hubble expansion of the Universe sug- 
gests that the density, velocity and gravitational fields may 
be decomposed into a sum of terms describing the uniform 
expansion and fluctuations against the background: 
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m p, u 



a(t) , . 



(A.ll) 



The term a[t)u is called the peculiar velocity. In cosmology, 
one also often employs the density contrast defined as (5 = 
p — 1, which gives the fluctuation against the normalized 
background density. Taking p, it, and </3g as functions of the 
comoving coordinate x = r/a{t) and using (lA.SL IIA.6I 1 and 
l|A.7fl . we rewrite the Euler-Poisson system in the form 



a 1„ ~ 

-2—u Vx^g, 

a a 



0, 

AnGgo 



ip- 



(A.12) 
(A.13) 
(A. 14) 



dtu + [u ■ Vx)u 
dtp + \7x ■ (pu) 

"x'fig 

Note the Hubble drag term —2{a/a)u in the right-hand side 
of IIA.I2II representing the relative slowdown of peculiar ve- 
locities due to the uniform expansion. 

around the trivial 



Formally linearizing IIA.12|l - l|A.14|l 

zero solution, one obtains the following ODE for the linear 
growth factor T(t) of density fluctuations: 

^(aV) =47tG^o^. (A. 15) 

The only solution of this equation that stays bounded (in- 
deed, vanishes) at small times is usually referred to as the 
growing mode. As we shall shortly see, it is convenient to 
choose the amplitude factor r of the growing mode to be 
a new 'time variable,' which in an Einstein-de Sitter uni- 
verse is proportional to t^^^ . It is normalized such that 
To = r(to) = 1- Rescaling the peculiar velocity and the grav- 
itational potential according to 



~ AtiGqot 
fs = 'Ps 



(A.16) 



and using the fact that in an Einstein-de Sitter universe 
dln(a^f)/dr = 3/(2r), we arrive at the following form of the 
Euler-Poisson system, which we use throughout this paper: 



drV + {v ■ \7x)v 

drP + Vx ■ {pv) 



--(. + v. 



0, 

p- 



(A.17) 
(A.18) 
(A.19) 



Suppose initially, i.e. at r = 0, a mass element is located 
at a point with the comoving coordinate q. Transported by 
the peculiar velocity fleld in the comoving coordinates, this 
element describes a trajectory x[q, r). Using the Lagrangian 
coordinate q to parametrize the whole continuum of mass 
elements, we recast 



A. 1711 and 1A.19II in the form 

f V^(^g) , (A.20) 
= i [(detVqO;)"^ - 1] . (A.21) 



DIx = -^{DrX 

It 



V72 



The density and peculiar velocity in Lagrangian variables 
are given by 



p(a;(q,r),T) = (det VqCc)" 
v{x{q,T),T) = DrX{q,T), 



(A.22) 



which automatically 
law llA.l 



satisfy the mass conservation 
Here Dt is the operator of Lagrangian 
time derivative, which in Lagrangian variables is the usual 
partial time derivative at constant q and in Eulerian 



variables coincides with the material derivative dr + v ■ Va, . 
The notation Va, in Lagrangian variables stands for the 
a;(q, T)-dependent differential operator with components 

= {dqj / dxi)\/ q ■ , which expresses the Eulerian gradient 
rewritten in Lagrangian coordinates, using the inverse 
Jacobian matrix. Note that Vx and Dt do not commute 
and that terms with Vx in the Lagrangian equations are 
implicitly non- linear. 

In one dimension, l|X 
quence: 



has an interesting conse- 



(A.23) 



Indeed, in one dimension IIA.211 takes the form 

V^^g = i[(V,x)-^-l]. (A.24) 

Multiplying this equation by VqX and expressing the flrst 
of the two ^-derivatives acting on <^g as a g-derivative, we 
obtain 



V, (Vx<(5g) = \7q 



(A.25) 



Eq. lA.23t is obtained from IIA.2511 by integrating in q. The 
absence of an arbitrary r-dependent constant is established 
either by assuming vanishing at large distances of both ipg 
and of the displacement x — q or, in the space-periodic case, 
by assuming the vanishing of period averages. 

1IA.23I1 to eliminate the (fie term in lfO(H and 

q, we 



Using ljA.23|l to eliminate the (fig term in 
introducing the notation ^ for the displacement x - 
obtain 



(A.26) 



The only solution to this equation that remains well-behaved 
for T ^ is the linear one ^ oc r. This solution has the two 
terms on the right-hand side of the one-dimensional version 
of ljA.20|l cancelling each other and hence gives a vanishing 
'acceleration' D^x. 

An approximate vanishing of acceleration takes place in 
higher dimensions as well. For early times, the Lagrangian 
map x{q,T) stays close to the identity, with displacements 
4(q,r) = x{q,T) — q small. Linearizing 1A.201 and 1A.211 
around zero displacement, we get the system 



^I'Pg 



(A.27) 
(A.28) 



Vq and detVqCC 



Here we use the fact that 
1 -I- Vq ■ ^. Using ljA.28|l to eliminate (fig in l|A.27^ . we get 
for 9 = \7q ■ £, an equation that coincides with ljA.26^ up 
to the change of variable ^ 1— > 6. Choosing the well-behaved 
linear solution for 9, solving for ^ and using the above ar- 
gument to eliminate a r-dependent constant, we see that, 
in the linearized equations, terms in the right-hand side of 
l|A.27fl cancel each other and the acceleration vanishes. This 
simplification justifies using the linear growth factor t as a 
time variable. 
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APPENDIX B: HISTORY OF MASS 
TRANSPORTATION 

The sub ject o f mass transportation was started by Gaspard 
Monge l|l781l l in a paper^'^ entitled Theorie des deblais et 
des remblais (Tiieory of cuts and fills) whose preamble is 
worth quoting entirely (our translation): 

When earth is to be moved from one place to another, the 
usage is to call cuts the volumes of earth to be transported and 
fills the space to be occupied after transportation. 

The cost of transporting one molecule being, all things other- 
wise equal, proportional to its weight and to the distance [espace] 
travelled and consequently the total cost being proportional to 
the sum of products of molecules each multiplied by the distance 
travelled, it follows that for given shapes and positions of the cuts 
and fills, it is not indifferent that any given molecule of the cuts 
be transported to this or that place in the fills, but there ought 
to be a certain distribution of molecules of the former into the 
latter, according to which the sum of these products will be the 
least possible, and the cost of transportation will be a minimum. 

Although clearly posed, the 'mass transportation prob- 
lem' was not solved , in m ore than one dimension, until 
Leonid Kantorovich il942l) formulated a 'relaxed' version, 
now called the Monge-Kantorovich problem: instead of a 
'distribution of molecules of the former mto the latter,' he 
allowed a distribution in the product space where more than 
one position in the fills could be associated with a position 
in the cuts and where the initial and final distributions are 
prescribed marginals (see Section 13.31 . In cosmospeak, he 
allowed multi-streaming with given initial and final mass 
distributions. Using the techniques of duality and of lin- 
ear programming that he had invented (see Appendix IC2II . 
Kantorovich was then able to solve the mass transporta- 
tion problem in this relaxed formulation. The techniques 
developed by Kantorovich found many applications, notably 
in economics, which in fact was his original motivation (he 
was awarded, together with T.C. Koopmans, the 1975 Nobel 
prize in this field). 

Before turning to more recent developments we must 
say a few words about the history of the Monge - Ampere 
equat ion. It was considered for the first time by lAmperel 
l)l82n'l for an unknown function z{x,y) of two scalar vari- 
ables. The equation is to be found on p. 65 of Ampere's huge 
(188 pages) mathematical memoir in the form 



Hr + 2Ks + Lt + M + N{rt - s^) = 0, 



(B.l) 



where in modern notation r — d^z/dx^, s — d^z/{dxdy), 
t = d^z/dy^ , and H, K, L, M, N are functions of x, y, z and 
the two first-order derivatives p = dz/dx and q = dz/dy. 
This extends the earlier work by lMong3 i 17841 see p. 126) 
concerning the equation without the Hessian term (TV = 0). 
Both Ampere and Monge were interested in methods of ex- 
plicit integration of these equations. Ampere also pointed 
out the way the equation changes under Legendre transfor- 
mations but there is no physical interpretation in terms of 
Lagrangian coordinates.^^ There is evidence that until the 



beginning of the 20th century the scientific community at- 
tributed the equation wit h the Hessian solely to Ampere 
ijBouil l)l86l p. 186) and IWebeil l|l90nl . p. 367)). But the 
joint attr ibution of |B.1|I to 'Monge and Ampere' is already 
found in llGoursad lT896). 

The subjects of mass transportation and of the Monge- 
Ampere equation came together when one of us (YB) showed 
the equivalence of the elliptic Monge-Ampere equation and 
of the mass transportation problem with quadratic cost: 
when initial and final distributions are non-singular, the 
optimal solution is actually one-to-one, so that nothing 
is los t by the Kantorovich relaxation trick (|Brenieij Il987l . 
I l99l|). For an extension o f this result to general costs see 
iGangbo fc McCannI il996t) ; a r eview of the man y recent pa- 
pers on the subject is given bv lAmbrosiol (|2r 



APPENDIX C: BASICS OF CONVEXITY AND 
DUALITY 

CI Convexity and the Legendre transformation 

A convex body may be defined by the condition that it co- 
incides with the intersection of all half-spaces containing it. 
Obviously, it is sufficient to take only those half-spaces lim- 
ited by planes that touch the body; such planes are called 
supporiing. 

Take now a convex function /(q), so that the set of 
points in the (3-1- l)-dimensional (q, /) space lying above its 
graph is convex. It follows that we can write 



f{q) = max xq- ,f*(x). 



(C.l) 



where the expression x ■ q — f* [x) specifies a supporting 
plane with the slope x for the set of points lying above the 
graph of / (see Fig. IGll for the one-dimensional case). The 
function f*{x), which specifies how high one should place a 
supporting plane to touch the graph, is called the Legendre 
transform of /(q).^* 

From Eq. ICH follows the inequality (known as the 
Young inequality) 



fil) + f* (x) ^ X ■ q for all x, q, 



(C.2) 



where both sides coincide if and only if the supporting plane 
with the slope x touches the graph of / at q. This fact, to- 
gether with the obvious symmetry of this inequality, implies 
that 



max X ■ q — f{q) 

9 



(C.3) 



Thus, the Legendre transform of a convex function is itself 
convex and the Legendre transform of the Legendre trans- 
form recovers the initial convex function. 

If however we apply IC.lll to a nonconvex function /, we 
obtain a convex function /* , whose own Legendre transform 
will give the convex hull of /, the largest convex function 
whose graph lies below that of /. 



The author's name appears in this paper as 'M. Monge,' where 
the 'M.' stands for 'Monsieur.' 

•^^ According to the biography of Ampere by L. Pearce Williams 
in the Dictionary of Scientific Biography, Ampere's paper was 
written - after he had switched from mathematics to chemistry 



and physics — with the purpose of facilitating his election to the 
Paris Academy of Science; one can then speculate that his men- 
tion of the Legendre transformation was influenced by Legendre's 
presence in this academy. ^^^^^^^^^ 
It was introduced in the one -dimensional ca se bv lMandelbroi j 
(1.939) and then generalized bv lFenchell ^94^). 
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f(q) 



f*(x) 



subgradieni 



Figure CI. A convex function f(q) and the geometrical con- 
struction of its Legendre transform f*(x). Also illustrated is the 
subgradient of f{q) at a non-smooth point. 



When / is both convex and difTerentiable, l)C.2^ be- 
comes an equality for x = Vq/(q). If /* is also difTeren- 
tiable, then one also has q — Vxf'ix). This is actually 
Legendre's original definition of the transformation, which 
is thus limited to smooth functions. Furthermore, if the orig- 
inal function is not convex and thus has the same gradient at 
separated locations, Legendre's purely local definition will 
give a multivalued Legendre transform. (In the context of 
the present paper this corresponds to multi-streaming.) 

Not all convex functions are differentiable (e.g. f{q) = 
|q|). But the Young inequality can be employed to define 
a useful generalization of the gradient: the subgradient of / 
at q is the set of all x for which the equality in holds 
(see Fig. ICH . If / is smooth at q, then Vqf{q) will be the 
only such point; otherwise, there will be a (convex) set of 
them. 

If a convex function has the same subgradient at more 
than one point, the function is said to lack strict convexity. 
In fact, strict convexity and smoothness are complementary: 
lack of one in a convex function implies lack of the other in 
the Legendre transform. 

Fo r further backgrou nd on convex analysis and geome- 
try, see lRockafellaJ (Il970l) . 



C2 Duality in optimization 

Suppose we want to minimize a convex function ^{q) subject 
to a set of linear constraints that may be written in matrix 
notation as Aq = b (vectors q satisfying this constraint are 
called admissible in optimization parlance). We now observe 
that 



inf $((/) = inf sup ^{q) 

Aq = b „ J, 



X ■ {Aq - b). 



(C.4) 



Indeed, should Aq not equal b, the sup operation in x will 
give infinity, so such q will not contribute to minimization. 
Here we use the inf /sup notation instead of min/max be- 
cause the extremal values may not be reached, e.g., when 
they are infinite. 

Using IC.IL we rewrite this in the form 



inf sup y ■ q — ^*{y) — x ■ [Aq — b) 
= inf sup {y — A'^x) ■ q — <I>* (y) + x 



(C.5) 



where ^*(y) is the Legendre transform of ^{q) and is 
the transpose of A. Taking inf in q first, we see that the 
expression in the right-hand side will be infinite unless y = 
A^x. We then obtain the optimization problem of finding 



sup X ■ b — ^"(A x) 



(C.6) 



which is called dual to the original one. Note that there are 
no constraints on the dual variable x: any value is admissi- 
ble. 

Denoting solutions of problems l)C.4^ and ljC.6|l by q* 

and a;*, we see that 



- ^*{A^x*) - x* 



■6 = 0, 



(C.7) 



because the optimal values of both problems are given 
by ijC.Sfl and thus coincide. Furthermore, for any admissi- 
ble q and x 



Ho) 



-^"{A^x) 



xb>0. 



(C.8) 



because the right-hand sides of ljC.4fl and ljC.6|l cannot pass 
beyond their optimal values. 

Moreover, let equality IIC.711 be satisfied for some admis- 
sible q* and x*; then such q* and x* must solve the problems 
ljC.4|l and l|C.6^ . Indeed, taking e.g. x* for x in 

(EHJ and 

using 1C.7L we see that for any other admissible q 



$(q*) < ^q) 



(C.9) 



i.e., that q* solves the original optimization problem ljC.4|l . 

Convex optimization problems with linear constraints 
considered in this section are called convex programs. Their 
close relatives are linear programs, namely optimization 
problems of the form 



inf c - q — inf sup c - q — x ■ {Aq — b), 

Aq = b,q^O X 



(C.IO) 



where notation q ^ means that all components of the 
vector q are nonnegative. Proceeding essentially as above 
with c ■ q instead of $(q), we observe that in order not 
to obtain infinity when minimizing in q in 1C.5L we have 
now to require that A^x ^ c (i.e. c — Al'^ x ^ 0). The dual 
problem thus takes the form 



sup x ■ b 



(C.ll) 



with an admissibility constraint on x. Instead of l)C.7^ 
and ijcTst we obtain 



6 = c- 



{A'^x* - c) 



q* =0 



and 



X ■ b ^ c ■ q or {Ax- 



q < 0, 



(C.12) 



(C.13) 



the latter inequality being automatically satisfied for any 
admissible x,q. Note that for linear programs, the fact that 
l|(TT2)l holds for some admissible q* ,x* also implies that q* 
and x* solve their respective optimization problems. 

For further background on optimiza tion and duahty, 
see, e.g.. |Papadimitriou fc Steiglitj il982h . 
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C3 Why the analogue computer of Section 14.21 
solves the assignment problem 

We suppose that the analogue computer described in Sec- 
tion 14.21 has settled into equilibrium, which minimizes its 
potential energy 



U 



N 

E 

i = l 



under the set of constraints 



(C.14) 



(C.15) 



for all Our goal is here is to show that the set of equi- 
librium forces fij, acting on studs between row and column 
rods, solves the original linear programming problem of min- 
imizing 



under constraints 

N N 



(C.16) 



(C.17) 



for all i,j and that in fact forces fij take only zero and 
unit values, thus providing the solution to the assignment 
problem. 

Note first that if a row rod Ai and a column rod Bj 
are not in contact at equilibrium, then the corresponding 
force vanishes [fij — 0); if they are, then fij 0. Take 
now a particular pair of rods Ai and Bj that are in contact. 
At equilibrium, the force fij must equal forces exerted on 
the corresponding stud by Ai and Bj. We claim that both 
these forces must be integer. To see this, let us compute 
the force exerted by Ai. This rod contributes its weight, 

possibly decreased by the force that it feels from other 
column rods that are in contact with Ai. Each of these takes 
— 1 (its 'buoyancy') out of the total force, but we may have 
to add the force it feels in turn from other row rods with 
which it might be in contact. Proceeding this way from one 
rod to another, we see that all contributions, positive or 
negative, are unity, so their sum fij must be integer. The 
same argument applies to rod Bj. 

Does this process indeed finish or, at some stage, do we 
come back at an already visited stud and thus end up in an 
infinite cycle? In fact, for general set of stud lengths C — Cij, 
the latter cannot happen, because otherwise an alternating 
sum of some subset of stud lengths would give exactly zero - 
a zero probability event for a set of arbitrary real numbers. 

Consider now a row rod Ai. It is in contact with one 
or more column rods, whose combined upward push must 
equilibrate the unit weight of Ai . Since any of the latter rods 
exerts a nonnegative integer force, it follows that exactly one 
of these forces is unity, and all the other ones are zero. A 
similar argument holds for any column rod Bj. 

We have thus shown that all fij in the equilibrium equal 
1 or 0. One can of course ignore the vanishing forces. Then 
each row rod Ai is supported by exactly one column rod Bj , 
and each Bj supports exactly one Ai. This defines a one- 
to-one pairing, and we are only left with a check that this 
pairing minimizes 1C.16II . 



Observe that pushing a column rod down by some dis- 
tance A and simultaneously increasing by A the length of 
all studs attached to this rod will have no effect on positions 
and constraints of all other rods, hence on the equilibrium 
network of contacts. Moreover, due to constraints 10.1711 . 
the corresponding change in coefficients Cij will not change 
the cost function JO. 1611 in any essential way, except of just 
subtracting A. 

We can use this observation to put all column rods at 
the same level, say at 2: = 0, adjusting dj to some new 
values c'ij. Thus, for every i, the row rod Ai rests on the stud 
with the largest height C — c'ij , so the equilibrium pairing 
maximizes the sum 

JV 

^(C-c^,)/,, (0.18) 
and thus minimizes 



APPENDIX D: DETAILS OF THE 
VARIATIONAL TECHNIQUE FOR THE 
EULER-POISSON SYSTEM 

In this appendix, we explain details of the variational pro- 
cedure outlined in Section |S| which proves that prescription 
of the density fields at terminal epochs r = and r = ro 
uniquely determines a regular and thus curl-free solution to 
the Euler-Poisson system lA.17t - l|A.19t . 



2 



The variational problem is posed for the functional 

TO 



dr 



3/2 



p\V\' 



(D.l) 



with four constraints: the Poisson equation 1A.19L which we 
repeat here for convenience, 

p-1 



(D.2) 



(D.3) 



the mass conservation 1A.18L also repeated here, 
drp + ■ (pv) = 0, 
and the two boundary conditions 

p(a;,0) = l and p{x,To) = po{x). (D.4) 

In the sequel, we shall always denote by J J the double in- 
tegration over ^ T ^ ro and over the whole space domain 
in X provided that the integrand vanishes at infinity suffi- 
ciently fast, or over the periodicity box in the case of periodic 
boundary conditions. A single integral sign J will always de- 
note the integration over the relevant space domain in x. 

First, we make this problem convex by rewriting the 
functional and constraints in a new set of variables with the 
mass fiux J{x,t) = p{x,t) v{x,t) instead of the velocity v. 
The mass conservation constraint, which was the only non- 
linear one in the old variables, becomes now linear: 



Those readers familiar with linear programming will recognize 
that the proof just presented is based on two ideas: (i) the total 
unimodularity of the matrix of constraints in terms of which the 
equalities in 1C.17I can be written and (ii) the com plementary 
slackness (see, e.g.. IPaoadimitriou SteiglitdllflSSl . sections 3.2 
and 13.2). 
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9rP + ■ J = 0, 



(D.5) 



and one can check that the density of kinetic energy takes 
the form 

p\v\ - 



2 
or 

2p 
where 



2p' 



max (pc + J ■ m) 

c, m: c+|mp/2^0 



= max (pc + J ■ m — F{c, m)), 



F{c,m) = I 



ifc+|m|V2<0 
+00 otherwise. 



(D.6) 



(D.7) 



Note that in l|D.6|l the variables c, m, as well as p, J, are 
functions of {x,t). The action functional may now be writ- 
ten as 



P 2 



3/2 ,3 



d X dr, 



(D.l 



and turns out to be convex. 

To see this, first note that the operation of integration 
is linear and thus preserves convexity of the integrand. The 
integrand is a positive quadratic function of Va:4> ^nd there- 
fore is convex in </!>; furthermore, (ID.6II implies that it is also 
convex in {p,J), since the kinetic energy density |Jp/2p is 
the Legendre transform of the function F{c, m), which itself 
is convex. 

Note also that by representing the kinetic energy den- 
sity in the form , we may safely allow p to take negative 
values: the right-hand side being in that case +oo, it will not 
contribute to minimizing ID.H . 

We now derive the dual optimization problem. We intro- 
duce the scalar Lagrange multipliers tp{x,t), i?in(a;), i?o(a;) 
and d{x,t) for the Poisson equation (ID.2II . the boundary 
conditions ljD.4^ . and the constraints of mass conservation 
, respectively, and observe that the variational problem 
may now be written in the form 



inf sup 

c+|Trep/2<0 



d a; dr 



f 3 
+Q(dTp + ■ J) + T^''^ \ pc^J-m + -j V^<i 



(D.9) 



/ i9i„(x)(p(a;,0) - l)d=' 



i9o(a;)(p(a;,ro) - po(a;))d x. 



To see that ID.Qt is indeed equivalent to minimizing lU.H 
under the constraints (ID.3II or ID.SL 1D.2II . and IID.411 . ob- 
serve that for those p, J,0 that do not satisfy the con- 
straints, the sup operation over d, ip, i^in, i^o will give pos- 
itive infinity; the sup will be finite (and thus contribute to 
the subsequent minimization) only if all constraints are sat- 
isfied. (This argument is the functional version of what is 
explained in Appendix IC2I for the finite-dimensional case.) 
Performing an integration by parts in the r variable in 
and using the boundary conditions on the mass density 
JD.4L we find that i?in(a;) = 9{x,0) and i?o(a;) = d{x,To). 
Integrating further by parts in the x variable, assuming that 
boundary terms at infinity vanish (or that we have periodic 
boundary conditions in space) and rearranging terms, we get 



inf sup 



• dr (c 



.3/2 



2r 



^) 



4t-3/2 



(D.IO) 



4^3/2 



(a;,0)d^a;+ / 9{x,to) pq{x)A'^x. 



Performing minimization with respect to p,J,4' first, as in 
l(a5|l of Appendix IC2l we see that the following two equali- 
ties must hold (remember that p need not be positive at this 
stage): 



™ = 7371^-^' 



(D.ll) 



so that terms linear in p and J vanish in ljD.10|l . It follows 
that c and m are determined by 9 and and that the con- 
straint c + |m|^/2 can be written 



(D.12) 



Also, the inf with respect to 4> is straightforward and gives 

r^'''V^(^g = VccV- (D-13) 

Using lITTTTl and iTTTsl in ilTToll . we arrive at the opti- 
mization problem of maximizing 



4r3/2 



+ 



9{x,To) po{x)d'^x - / e{x,0)d^x 



(D.14) 



under constraint llD.12ll . Eqs. <D.14t and IIP. 121 constitute 
a variational problem dual to the original one. 

As both the original and the dual variational problems 



have the same saddle-point formulation ljD.9|l or l|D.10^ , the 



optimal values of the two functionals l|D.l|l and ljD.14|l are 
equal. Let (p, J, ipg) be a solution to the original variational 
problem and 9,tp he a solution to the dual one. Subtracting 
the (equal) optimal values from each other, we may now 
write, similarly to IIC.711 . 



+ 



3/2 ■^t'^/^ 

^^|V.7M"-|-7/-)d^a;dr 



(D.15) 



+ I 9{x,0)d^x- I e{x,To)po{x)d^x^O. 



We are going to show that the left-hand side of lID.lSIl may 
be given the form of a sum of three nonnegative terms, each 
of which will therefore have to vanish. First, we rewrite the 
last two integrals, using the mass conservation constraint 
l|D.5^ and integrations by parts, in the form 



ar(6lp)d='a;dr = - / I {drO p + ^ ^9 ■ J) d?x dr. 



Second, we note that 



© 0000 RAS, MNRAS 000, 000-000 



Reconstruction of the early Universe as a convex optimization problem 25 



which follows from the Poisson constraint 1D.2II . Taking all 



this into account in l|D.15fl . we get, after a rearrangement of 
terms, 

^3/2 2 

J - V, 

p 



2r3/2 

+ / [ -p[dre + 



_i_,V.„2 + |^)d='xdr (D.16) 



4-^3/2 



3/2v7 



\/^ip)\'^ d^xdr = 0. 



The left-hand side is a sum of three nonnegative terms (the 
second is so by IIP. 12^ '). all of which must thus vanish. This 
gives 

1 _ „ _ 1 



P 



r3/2 



r3/2 



and 



dre + 



2t3/^ 



2r^ 



0, 



(D.17) 



(D.18) 



wherever p is non-vanishing (otherwise the left-hand-side is 
non-positive by lD.12t 'l. The last equality turns into the 
Euler equation 

drV + {v-\/^)v = -^{v + \/^iPs) (D.19) 
At 

by taking the gradient and using lD.17t . 

By l|D.17^ and l|D.18fl . any two hypothetically different 
minimizing solutions for either variational problem give rise 
to the same velocity potential and to the same gravitational 
potential (up to insignificant constants) and thus define the 
same solution (p, v,ifg) to the Euler-Poisson equations with 
the boundary conditions iPAl and the condition of curl-free 
velocity. 

Moreover, for any such solution {p,v,ifg), one can 
use (ID. 1711 to define 9 and tp that satisfy ID. 1811 and 
thus 



(|DTT2j. By ljD.16^ . the values of functionals / and I 
evaluated at these functions will coincide; together with con- 
vexity this implies, by an argument similar to that given 
in Appendix IC2I concerning ljC.9|l . that such {p,v,(pg) and 
{9, ip) in fact minimize both functionals under the corre- 
sponding constraints. 

This means that a (curl-free) velocity field, a gravita- 
tional field and a density fields (D,(/3g,p) will satisfy the 
Euler-Poisson equations (IA.17t - l|A.19t (repeated as (ID. 1911 . 
ljD.3^ . and ljD.2^ in this Appendix) and the boundary con- 
ditions (|D.4fl if and only if they minimize (ID.lfl under the 
corresponding constraints. This establishes uniqueness. 
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